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Summary 

We  report  here  on  the  work  performed  during  the  third  year  (october  2011  - 
October  2012)  of  the  contract  FA  8655-10-C-4002  on  Multiscale  problems  in 
materials  science:  a  mathematical  approach  to  the  role  of  uncertainty. 

We  recall  that  the  bottom  line  of  our  work  is  to  develop  affordable  numer¬ 
ical  methods  in  the  context  of  heterogeneous,  possibly  stochastic  materials. 
Many  partial  differential  equations  of  materials  science  indeed  involve  highly 
oscillatory  coefficients  and  thus  small  length-scales.  When  the  microstructure 
of  the  materials  is  periodic,  or  random  and  statistically  homogeneous,  homog¬ 
enization  theory  can  be  used,  and  allows  to  appropriately  define  averaged 
equations  from  the  original  oscillatory  equations.  The  theoretical  aspects  of 
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these  problems  are  now  well-understood,  in  general.  On  the  other  hand,  the 
numerical  aspects  have  received  less  attention  from  the  mathematics  commu¬ 
nity,  in  particular  in  the  case  of  stationary  ergodic  random  problems,  which 
are  one  instance  often  used  for  modelling  uncertainty  of  continuous  media. 
In  that  latter  case,  standard  methods  available  in  the  literature  often  lead  to 
very,  and  sometimes  prohibitively,  costly  computations. 

The  situation  is  even  more  challenging  when  no  structural  assumption 
(periodicity,  statistical  homogeneity,  . . . )  on  the  materials  microstructure 
can  be  made.  In  the  absence  of  such  an  assumption,  homogenization  theory 
still  holds,  but  does  not  provide  any  explicit  formulae  amenable  (even  pos¬ 
sibly  after  some  approximation)  to  numerical  computation.  One  possibility 
is  then  to  directly  address  the  original  problem  (rather  than  passing  to  the 
limit  of  inhnite  scale  separation),  and  to  use  dedicated  numerical  approaches 
for  such  multiscale  problems,  such  as,  for  instance,  the  Multiscale  Finite 
Element  Method  (MsFEM). 

In  this  report,  we  hrst  consider  a  variant  of  stochastic  homogenization, 
well  suited  to  model  materials  that  are  periodic  up  to  a  random  deformation. 
We  have  already  considered  this  variant  in  our  previous  report  [3],  but  with  a 
perspective  different  from  the  one  here.  This  variant  admits  a  homogenized 
limit.  However,  the  homogenized  matrix  is  expensive  to  compute,  as  is  often 
the  case  in  stochastic  homogenization.  We  propose  here  an  efficient  MsFEM 
type  approach  dedicated  to  that  setting. 

We  next  turn  to  studying  the  robustness  of  the  MsFEM  approach  to 
perturbations  of  the  equation  coefficients  that  are  non  oscillatory.  Our  idea 
is  that  MsEEM  approaches  are  devoted  to  capturing  the  highly  oscillatory 
modes  of  the  solution,  which  are  poorly  captured  by  a  standard  FEM  ap¬ 
proach  using  a  limited  number  of  degrees  of  freedom.  When  the  coefficient  in 
the  equation  is  modihed  by  a  non-oscillatory  component,  the  high  frequencies 
are  not  modified,  and  the  MsFEM  approach  can  be  expected  to  be  robust 
with  respect  to  these  perturbations.  This  is  exactly  the  question  we  consider 
in  the  second  part  of  this  report. 

The  works  described  below  have  been  performed  by  Claude  Le  Bris  (PI), 
Frederic  Legoll  (Co-PI)  and  Florian  Thomines  (third  year  Ph.D.  student). 
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1  Introduction 


During  this  third  year  of  contract,  we  have  pursued  our  effort  on  developing 
affordable  numerical  methods  in  the  context  of  stochastic  homogenization. 

Many  partial  differential  equations  of  materials  science  indeed  involve 
highly  oscillatory  coefficients  and  small  length-scales.  Homogenization  the¬ 
ory  is  concerned  with  the  derivation  of  averaged  equations  from  the  original 
oscillatory  equations,  and  their  treatment  by  adequate  numerical  approaches. 
Stationary  ergodic  random  problems  are  one  of  the  most  famous  instances  of 
mathematical  uncertainty  of  continuous  media. 

The  purpose  of  this  report  is  to  present  the  recent  progress  we  have  made 
on  this  topic,  with  the  aim  to  make  numerical  random  homogenization  more 
practical.  As  already  mentioned  in  our  two  previous  reports,  because  we 
cannot  embrace  all  difficulties  at  once,  the  case  under  consideration  here 
is  a  simple,  linear,  scalar  second  order  elliptic  partial  differential  equation 
in  divergence  form,  for  which  a  sound  theoretical  groundwork  exists.  We 
focus  here  on  the  different  manners  the  problem  can  be  handled  from  the 
computational  viewpoint. 

In  this  report,  we  are  concerned  with  various  questions  related  to  the 
MsFEM  approach.  We  recall  that  this  is  one  approach  (among  others,  see 
e.g.  [13]  for  an  alternative)  to  address  highly  oscillatory  problems  when  the  as¬ 
sumptions  needed  by  the  homogenization  theory  on  the  materials  microstruc¬ 
ture  (such  as  periodicity,  statistical  homogeneity,  . . . )  are  not  met.  We  have 
already  contributions  extending  the  range  of  applicability  of  that  approach, 
see  our  publication  [4]  and  the  previous  reports  [2,  Section  4]  and  [3,  Section 
3], 

We  begin,  in  Section  2,  with  a  brief  description  of  periodic  homogenization 
and  the  MsFEM  approach  in  a  deterministic  setting.  The  only  purpose  of 
that  section  is  the  consistency  of  this  report. 

In  Section  3,  we  consider  a  variant  of  stochastic  homogenization,  intro¬ 
duced  by  the  PI  and  co-workers  in  [9,  10]  some  years  ago.  This  model  is 
adequate  to  represent  materials  that  are  random  deformations  of  a  perfect 
periodic  material.  A  typical  example  is  a  composite  material  with  hbers. 
Fibers  are  all  identical,  they  would  be  located  on  a  periodic  lattice  in  the 
ideal  situation.  However  (for  instance  as  a  consequence  of  the  manufacturing 
process),  their  actual  positions  are  now  random.  In  our  previous  report  [3, 
Section  5],  we  have  presented  and  analyzed  a  procedure  to  practically  ap- 
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proximate  the  homogenized  matrix.  In  this  report,  we  propose  a  MsFEM 
type  approach  to  compute  an  approximation  of  the  solution  to  the  original 
highly  oscillatory  problem.  Although  the  setting  is  stochastic,  it  turns  out 
that  the  approach  we  propose  does  not  require  the  recomputation  of  the  Ms¬ 
FEM  basis  functions  for  each  new  realization  of  the  material.  This  is  why 
the  proposed  approach  is  much  less  expensive  than  the  natural  adaptation  of 
the  MsFEM  approach  to  the  stochastic  problem  at  hand.  The  performance 
of  the  approach  is  illustrated  by  some  numerical  tests,  which  demonstrate 
the  accuracy  of  the  computed  approximation. 

In  Section  4,  we  next  turn  to  a  question  which  is  originally  motivated  by 
inverse  problems  in  multiscale  science.  Assume  that  we  model  our  hetero¬ 
geneous  materials  with  the  oscillatory  coefficient  b{x)A^{x),  and  that  we  do 
not  have  a  good  knowledge  of  b.  This  situation  corresponds  to  the  case  when 
we  accurately  know  the  properties  of  our  materials,  up  to  a  slowly  varying, 
macroscopic  envelop.  Otherwise  stated,  the  high  frequency  modes  are  well 
identihed,  but  the  way  they  change  over  macroscopic  distances  is  not  well 
characterized.  One  possibility  to  dehne  a  relevant  b  is  to  search  for  b  so 
that  the  homogenized  properties  of  the  materials  are  as  close  as  possible  to 
the  ones  that  are  observed  in  practice.  We  thus  want  to  optimize  on  b,  and 
are  thus  going  to  iterate  on  this  function,  until  we  hud  the  best  one.  For 
each  iterate,  we  need  to  very  efficiently  compute  the  corresponding  homog¬ 
enized  properties  or  the  homogenized  solution  (because  such  a  computation 
is  needed  at  each  iteration  of  the  optimization  loop).  In  the  specihc  setting 
considered  here,  the  high  frequency  modes  are  independent  of  b,  since  b  is 
only  slowly  varying.  We  thus  expect  that  the  MsFEM  approach,  which  aims 
at  properly  taking  into  account  the  high  frequencies  present  in  the  problem, 
should  be  insensitive  to  the  choice  of  b.  Otherwise  formulated,  we  expect  the 
MsFEM  basis  functions  to  be  robust  with  respect  to  modihcations  on  b.  In 
Section  4.1,  we  consider  the  above  case  where  the  coefficient  reads  b{x)A^{x), 
and  indeed  show  that  the  MsFEM  basis  functions  can  be  computed  indepen¬ 
dently  of  b.  As  shown  in  Section  4.2,  the  situation  is  different,  and  more 
challenging,  when  the  coefficient  reads  b{x)  -|-  A^{x).  For  that  latter  case, 
we  propose  an  approximation  strategy  based  on  Proper  Generalized  Decom¬ 
position  ideas,  and  numerically  demonstrate  its  efficiency.  For  the  sake  of 
simplicity,  and  because  again  we  cannot  embrace  all  difficulties  at  once,  we 
only  consider  deterministic  models  in  that  Section  4. 
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We  eventually  collect  in  Section  5  some  possible  future  directions  of  re¬ 
search  that  we  may  consider  if  BOARD  decides  to  renew  our  funding. 

2  Periodic  homogenization  and  MsFEM  ap¬ 
proaches 

[Detailed  presentation  can  he  read  in  [1,  2].] 

For  the  consistency  of  this  report,  we  recall  in  this  section  some  ground¬ 
work  on  periodic  homogenization  and  on  related  approaches  for  deterministic, 
non  necessarily  periodic,  heterogeneous  materials.  More  details  can  be  read 
in  our  hrst  report  [2]  and  references  therein,  and  also  in  the  review  article  [1] 
that  we  published. 

The  problem  under  consideration  writes 

— div  [A^{x)Vu'^]  =  f{x)  in  V,  u^{x)  =  0  on  dV,  (1) 

where  "D  is  a  regular,  bounded  domain  in  and  where,  for  any  e,  the  matrix 
is  symmetric,  bounded  and  dehnite  positive.  The  parameter  e  encodes  the 
typical  size  of  the  heterogeneities.  We  manipulate  for  simplicity  symmetric 
matrices,  but  the  discussion  carries  over  to  non  symmetric  matrices  up  to 
slight  modihcations. 

2.1  Periodic  homogenization 

Assume  that,  in  (1),  the  matrix  A^  reads 

A^(a;)  =  Aper  (^)  (2) 

where  the  matrix  Aper  is  symmetric  dehnite  positive  and  ^Aperiodic.  In 
this  framework,  it  is  well  known  that,  as  e  — >  0,  the  solution  to  (l)-(2) 
converges  to  u*  solution  to  the  homogenized  problem 

—div  [Apg,.VM*]  =  f{x)  in  D,  u*{x)  =  0  on  dV,  (3) 

where  the  homogenized  matrix  A*  j,  reads 
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where  Q  =  (0, 1)'^  and  where,  for  any  p  G  the  so-called  corrector  Wp  is 
the  (unique  up  to  the  addition  of  a  constant)  solution  to 

-div  [Aper  {p  +  Vwp)]  =  0  on  M"*, 

Wp  is  Z'^-periodic.  ^ 

The  practical  interest  of  the  approach  is  evident.  No  small  scale  e  is  present 
in  the  homogenized  problem  (3).  At  the  price  of  only  computing  d  peri¬ 
odic  problems  (5)  (as  many  problems  as  dimensions  in  the  ambient  space), 
the  solution  to  (l)-(2)  can  be  efficiently  approached  for  e  small.  In  con¬ 
trast,  a  direct  attack  of  (l)-(2)  would  require  taking  a  meshsize  smaller  than 
e,  to  appropriately  capture  the  variation  of  the  materials  properties  at  the 
microstructure  scale.  The  difficulty  has  been  circumvented. 


2.2  MsFEM  approach 

The  homogenization  result  recalled  in  Section  2.1  heavily  relies  on  the  pe¬ 
riodicity  assumption  (2).  Although  it  is  possible  to  somewhat  relax  this 
assumption  and  still  obtain  explicit  formulae  for  the  homogenized  matrix, 
there  are  many  cases  of  practical  interest  for  which  the  existence  of  a  homog¬ 
enized  matrix  is  known,  but  no  explicit  formulae  are  available.  For  practical 
purposes,  alternative  approaches  are  needed. 

The  Multiscale  Finite  Element  Method  (MsFEM  approach)  is  one  such 
approach  (note  that  other  approaches  have  also  been  proposed  within  the 
same  paradigm,  we  refer  e.g.  to  [13]).  The  MsEEM  is  designed  to  directly 
address  the  original  problem  (1)  by  performing  a  variational  approximation 
using  pre-computed  basis  functions  yf  that  are  adapted  to  the  problem.  The 
method  is  not  restricted  to  the  periodic  setting,  in  contrast  to  the  homog¬ 
enization  theory  recalled  above.  We  do  not  assume  that  (2)  holds.  In  the 
sequel,  we  briefly  describe  the  approach,  and  refer  to  our  publication  [4]  (see 
also  [5])  for  more  details  and  comprehensive  numerical  tests. 

In  the  sequel,  we  argue  on  the  variational  formulation  of  (1): 

Eind  G  Hq{V)  such  that,  Vn  G  Hq{V),  =  b{v),  (6) 

where 

Ae{u,v)=  [  {Vv{x)fA^{x)Vu{x)dx  and  b{v)  =  [  f{x)v{x)dx. 

Jv  Jv 

We  introduce  a  classical  PI  discretization  of  the  domain  V,  with  L  nodes, 
and  denote  y^,  i  =  1, ...  ,L,  the  basis  functions. 
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Definition  of  the  MsFEM  basis  functions  Several  definitions  of  the 
MsFEM  basis  functions  have  been  proposed  in  the  literature  (see  e.g.  [16, 
14,  11]).  They  give  rise  to  different  variants  of  the  method.  In  the  following, 
we  present  one  of  these  variants.  For  any  finite  element  (e.g.  triangle)  K,  we 
consider  the  problem 

-div  =0  in  K, 

=X°Ik  ondK. 

By  construction,  these  functions  that  are  numerically  precomputed, 

encode  the  fast  oscillations  present  in  (1). 

Note  the  similarity  between  (7)  and  the  corrector  problem  (5).  Note  also 
that  the  problems  (7),  indexed  by  K,  are  all  independent  from  one  another. 
They  can  hence  be  solved  in  parallel,  using  a  discretization  adapted  to  the 
small  scale  e. 

Macro  scale  problem  We  now  introduce  the  finite  dimensional  space 

Wh  :=  span{yf,  i  =  1,...,L}, 

where  yf  is  such  that  yflj^  =  all  K,  and  proceed  with  a  standard 

Galerkin  approximation  of  (6)  using  Wft: 

Find  G  Wh  such  that,  Vn  G  Wh,  Ae{vf,,v)  =  b{v).  (8) 

The  function  is  the  MsFEM  approximation  of  the  exact  solution  .  Note 
that  the  dimension  of  Wh  is  equal  to  L\  the  formulation  (8)  hence  requires 
solving  a  linear  system  with  only  a  limited  number  of  degrees  of  freedom. 

Numerical  illustration  In  order  to  illustrate  the  MsFEM  approach,  we 
solve  (1)  in  a  one  dimensional  setting  with 

A^{x)  =  5  +  SOsin^  ^ 

on  the  domain  V  =  (0, 1),  with  e  =  0.025  and  /  =  1000.  We  subdivide  the 
interval  (0, 1)  in  L  =  10  elements.  On  Figure  1,  we  plot  the  MsFEM  basis 
functions  in  a  reference  element  and  the  MsFEM  solution 
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Figure  1:  Left:  Multiscale  basis  functions  fhe  reference  element. 

Right:  MsFEM  solution  in  the  domain  (0, 1). 


3  A  MsFEM  type  approach  for  a  variant  of 
stochastic  homogenization 

[Work  expanded  in  [5].] 

As  announced  in  the  conclusion  of  our  previous  report  [3],  we  consider 
here  a  variant  of  the  classical  setting  of  stochastic  homogenization,  originally 
introduced  a  few  years  ago  in  [9,  10],  and  propose  for  that  variant  a  MsFEM 
type  approach.  We  briefly  review  the  problem  under  consideration,  before 
recalling  the  corresponding  homogenization  results  and  eventually  describing 
our  contribution. 


3.1  A  variant  of  stochastic  homogenization  and  its  ho¬ 
mogenized  limit 

The  equation  under  consideration  is 


-div 


X 


=  f{x)  in  a;)=0  on  dV,  (9) 


where  the  matrix  A  is  the  composition  of  a  Z'^  periodic  matrix  Aper  with  a 
stochastic  diffeomorphism  <F: 


A 


:=  A 


per 


(10) 


We  assume  that,  almost  surely,  the  map  <h(-,a;)  is  a  well-behaved  diffeomor¬ 
phism  from  to  (in  the  sense  that  Esslnf^^f^  (det(V<F(a:,  a;)))  =  z/  > 
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0  and  EssSup^gQ ci;)|  =  M  <  +cxd),  and  that  it  satisfies 

V<h  is  stationary.  (11) 

Formally,  such  a  setting  is  well  suited  to  model  materials  that  are  periodic, 
in  a  given  reference  configuration.  The  latter  is  only  known  up  to  a  cer¬ 
tain  randomness.  Materials  we  have  in  mind  are  ideally  periodic  materials, 
where  some  random  deformation  (modelled  by  <l>)  has  been  introduced,  for 
instance  during  the  manufacturing  process.  Assumption  (11)  means  that  V<h 
is  statistically  homogeneous,  i.e.  the  randomness  is  the  same  anywhere  in 
the  material.  We  refer  e.g.  to  [2,  Section  2.2]  for  a  brief  discussion  of  the 
notion  of  stationarity  in  the  context  of  random  homogenization. 

The  problem  (9)-(10)  admits  a  homogenized  limit  when  e  vanishes.  It  is 
indeed  shown  in  [9]  that,  under  the  above  assumptions,  the  solution 
to  (9)-(10)  converges  as  e  goes  to  0  to  u*,  solution  to  the  deterministic  ho¬ 
mogenized  problem 

— div  =  f{x)  in  V,  u*{x)  =  0  on  dV. 

The  homogenized  matrix  A*  is  given  by,  for  any  1  <  i,  j  <  d. 


A*  =  det  (  E  I  /  V<I>(|/,  ■)dy 


'Q 


-1 


X 


E 


ejAper  ($  ^  {y,  ■))  {ej  +  Vwejiy,  ■))  dy  )  , 


where  Q  =  (0, 1)^  and  where,  for  any  p  G  Wp  solves  the  corrector  problem 


"  -div  [Aper  (<h  ^{y,uj))  {p  +  Vwp{y,u))]  =  0  in  R'^, 
Wp{y,u)  =  Wp  (^^~^{y,uj),uj)  ,  Vwp  is  stationary, 

e(  [  Vwp{y,-)dy)  =0. 

\d$(Q,.)  / 


(12) 


3.2  A  MsFEM-type  approach 

Our  aim  here  is  to  propose  a  MsFEM-type  approach  for  (9)-(10).  One  mo¬ 
tivation  is  that,  as  is  standard  in  stochastic  homogenization,  the  corrector 
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problem  (12)  is  set  on  the  complete  space  and  is  thus  challenging  to  solve 
in  practice. 

Note  that  efficient  MsFEM  approaches  are  not  easy  to  derive  in  stochastic 
settings,  when  the  equation  of  interest  writes  in  the  general  form 


—div[A^{x,u)Vu'^{x,uj)]  =  f{x)  inV,  u^{x,u)  =  0  ondV.  (13) 


As  pointed  out  in  our  Erst  report  (see  [2,  Section  4.1]),  a  natural  adapta¬ 
tion  of  the  deterministic  MsFEM  approach  presented  in  Section  2.2  to  the 
problem  (13)  would  involve  computing,  for  each  new  random  realization  of 
the  matrix  A^{x,uj),  new  highly  oscillatory  basis  functions  (see  indeed  (7)). 
This  is  prohibitively  expensive.  However,  in  particular  settings,  dedicated 
MsFEM-type  approaches  can  be  proposed.  We  have  considered  in  our  pre¬ 
vious  reports  (see  [2,  Section  4.2]  and  [3,  Section  3])  a  weakly-stochastic 
setting,  where  the  random  matrix  A^{x,u)  in  (13)  is  a  small  perturbation  of 
a  deterministic  matrix,  and  proposed  for  that  particular  setting  an  appropri¬ 
ate  and  efficient  MsEEM  approach  (see  our  publication  [4]).  In  what  follows, 
we  consider  the  particular  setting  (9)-(10),  and  use  in  an  essential  manner 
the  fact  that  it  is  built  upon  a  periodic  matrix,  randomly  deformed. 

We  know  from  (12)  that  the  expectation  of  Wp  is  a  periodic  function. 
Our  approach  is  based  on  approximating  the  corrector  Wp  in  (12)  by  a  periodic 
function  fcP®''. 

To  proceed,  it  is  useful  to  write  the  corrector  problem  (12)  in  a  variation- 
nal  form.  As  shown  in  [9],  we  have  that 


E 


,■) 


{V'4){y,u)Y Ap^,  ($  ^{y,uj))  {p +  Vwp{y,u))dy 


0 


for  all  Ip  stationary,  and  where  ip  =  ip  o  ^  The  above  expression  can  be 
rewritten,  after  a  change  of  variables,  as 


E 


det(V4>) 


Uq 


Vipy  A. 


per 


(p-h(V<h)  '^Vwp) 


=  0. 


We  introduce  the  notation  <F  =  E(<l>).  Our  idea  is  to  approximate  the  random 
function  Wp  by  a  periodic  function  solution  to 


Aper(p+(V<h)“^VhiP“)  =0 


(14) 
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for  all  functions  tjj  that  are  Z'^-per iodic.  Note  that  is  uniquely  dehned 
(up  to  the  addition  of  a  constant)  by  the  above  problem. 

Remark  1  In  general,  the  funetion  Wp  is  not  periodic.  An  explicit  counter¬ 
example  is  given  in  [9],  In  the  numerical  tests  below,  we  will  consider  that 
particular  example,  and  show  that  our  approach  yields  accurate  results  even 
in  that  dijficult  case.  See  also  Remark  4  below. 

Definition  of  the  basis  functions  As  in  Section  2.2,  we  introduce  a 
classical  PI  discretization  of  the  domain  T>,  with  L  nodes,  and  denote 
i  =  1, . . . ,  L,  the  basis  functions.  We  denote 

Vh  :=  Span(x°) 

the  associated  hnite  dimensional  space. 

Let  tyP®''  be  a  solution  to  (14).  We  set 

:=  uiP®"  (4>"^(a;,o;))  (15) 

and  introduce  the  vector  W {x,  u)  G  whose  components  are  given  by 

Wj{x,u)  =  eJW{x,uj)  :=  t(;gPP(a;,a;),  1  <  j  <  d. 

The  highly  oscillatory  basis  functions  are  dehned  by 

Xt{x,uj)  :=Xi  [x  +  eW  ,  1  <  i  <  L. 

Remark  2  In  the  case  when  <h(a;,  cn)  =  x,  the  problem  (9)-(10)  under  con¬ 
sideration  is  exactly  the  highly  oscillatory  problem  (l)-(2),  with  a  periodic 
matrix  coefficient.  In  that  case,  the  approach  proposed  above  is  identical  to 
the  MsFEM-type  approach  proposed  in  [11]  to  address  (l)-(2). 

Macro  scale  problem  We  introduce  the  hnite  dimensional  space 

Wh  ■=  Span(xf) 

and  proceed  with  a  Petrov-Galerkin  approximation  of  (9)-(10).  The  numer¬ 
ical  approximation  u\  G  Wh  is  dehned  as  the  unique  solution  to  the  weak 
formulation 

Vn  G  Vh,  j  {yv{x)Y  [^,uyjVul{x,u)dx  =  j  f{x)v{x)dx. 

(16) 
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The  main  feature  of  the  proposed  approach  is  that,  although  the  highly 
oscillatory  basis  functions  yf  are  stochastic  (and  thus  depend  on  the  realiza¬ 
tion  of  the  random  material),  we  actually  do  not  have  to  solve  a  new  problem 
for  each  new  realization  of  the  material  (i.e.,  for  each  new  realization  of  the 
diffeomorphism  $)  to  compute  these  basis  functions.  The  basis  functions  are 
indeed  given  by  (15),  where  the  deterministic  function  has  been  pre¬ 
computed.  For  each  new  realization,  we  thus  only  have  to  evaluate  the  new 
basis  functions.  This  is  the  main  advantage  in  terms  of  cost  in  comparison 
with  a  natural  application  of  the  MsFEM  approach  on  the  problem  (9)-(10). 

Note  that,  for  each  new  realization  of  the  material,  we  have  to  recompute 
the  stiffness  matrix  of  the  problem,  which  is  given  by 

=  j  (Vx°(a;))^Aper  '^X%x,u)dx,  1  <  i,j  <  L. 

A  natural  adaptation  of  the  MsFEM  approach  on  the  problem  (9)-(10)  would 
also  involve  recomputing  the  stiffness  matrix  for  each  new  realization. 

Remark  3  We  have  chosen  in  (16)  to  perform  a  Petrov- Galerkin  approxi¬ 
mation  of  the  problem,  where  the  space  Vh  of  the  test  functions  v  is  different 
from  the  space  Wh  of  the  numerical  solution  uf.  It  is  also  possible  to  use  a 
Galerkin  approximation,  which  yields  results  the  accuracy  of  which  is  similar, 
although  not  as  good,  as  the  results  presented  below. 


Some  elements  of  analysis  in  a  one-dimensional  setting  In  the  one¬ 
dimensional  setting,  the  corrector  problem  (12)  reads 


A. 

dy 


dlper  ($  ^(y,^))  (  1  + 

dw 


=  0  in 


w{y,uj)=w(^^  ^(y, cu), cu)  ,  is  stationary, 

eI/  >.-M!/)=0. 

dy  J 


This  problem  can  be  analytically  solved,  and  we  obtain  that 

-  1, 


dw  C 


(17) 
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where  C*  is  a  deterministic  constant  given  by 


1 

C 


E 


E  $'(|/,  Ody)  Aer(y) 


Since  w{y,uj)  =  w{^{y,u),u),  we  compute  by  the  chain  rule  that 


dw 

dy 


dw 

dy 


{(l>{y,uj),uj) 


d^ 


(18) 


(19) 


We  have  now  completely  characterized  the  solution  to  the  corrector  prob¬ 
lem  (17).  We  next  turn  to  the  problem  (14)  that  we  introduced  in  our 
MsFEM-type  approach.  In  the  one-dimensional  setting,  this  problem  reads 


dip  (  f  d^\  ^  dw'^‘^'^\ 

'o  ^  ^  \^)  ~d^ ) 


=  0 


for  all  functions  ip  that  are  E-periodic.  Again  using  the  specihcities  of  the 
one-dimensional  setting,  we  can  analytically  solve  this  problem,  and  obtain 

that  _ 

dfcP®''  d^  P  C 

dy  dy 

where  the  constant  C  is  again  given  by  (18).  Comparing  (19)  and  (20),  we 
see  that  the  exact  corrector  w  and  our  approximate  solution  ■uIp®'’  are  related 
by 

^~per  ^  f 

dy  \dy  J  ' 

This  somehow  shows  the  consistency  (at  least  in  the  one-dimensional  setting) 
of  our  approximate  problem  (14).  Note  however  that  the  one-dimensional 
case  may  be  misleading  in  that  respect,  as  it  is  the  only  case  where  the 
gradient  of  the  corrector  Wp  solution  to  (12)  is  of  the  form  of  “  a  periodic 
function  composed  with  the  diffeomorphism  In  general,  this  is  not 

the  case,  as  explicitly  pointed  out  in  [9].  Numerical  tests  are  therefore  of 
paramount  importance  to  validate  the  approach. 


Numerical  tests  We  have  considered  the  problem  (9)-(10)  on  the  domain 
V  =  (0, 1)^,  for  two  test  cases  represented  on  Figure  2.  In  both  cases,  the 
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periodic  matrix  Ap^i  represents  hard  inclusions  in  a  soft  material.  In  the 
hrst  test  case  (top  row  of  Figure  2),  inclusions  have  a  circular  shape,  and 
the  diffeomorphism  corresponds  to  a  translation  of  these  inclusions  by 
a  random  vector.  More  precisely,  on  each  cell  k  +  Q  (with  Q  =  (0, 1)'^),  the 
function  is  a  translation  by  a  random  vector  Xk{uj)  G  The  random 
variables  Xk  are  independent  and  identically  distributed.  In  the  second  case 
(bottom  row  of  Figure  2),  inclusions  have  a  T  shape,  and  the  diffeomorphism 
corresponds  to  a  rotation  of  these  inclusions  by  a  random  angle  0,  which 
can  take  four  values,  with  equal  probability:  6*  =  0,  7r/2,  vr  or  37r/2. 

Remark  4  The  second  test  case  is  inspired  by  (and  the  resulting  function  <I> 
is  very  close  to)  the  counter-example  discussed  in  [9],  It  is  shown  there  that, 
for  this  counter-example  case,  the  gradient  of  the  corrector  Wp  solution  to  (12) 
is  not  of  the  form  of  “  a  periodic  function  composed  with  the  diffeomorphism 
4)“^  ”.  Despite  this  fact,  we  show  below  that  the  ansatz  (15)  (and  the  resulting 
MsFEM-type  approach  described  above)  actually  yields  accurate  results  (see 
the  second  line  of  Table  1 ). 


Figure  2:  Left:  the  periodic  material  modelled  by  Right:  a  realization 
of  y4per(*F~^).  Top  row:  $  is  a  translation  of  circular  inclusions.  Bottom  row: 
<h  is  a  rotation  of  T-shaped  inclusions. 
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We  work  with  the  parameters  e  =  0.025  and  H 
between  two  random  fnnctions  ui  and  U2  is  dehned  by 


e{ui,U2)  =  E 


/  II  Vmi  —  Vm2||l2(x>)  \ 

V  ||Vn2||L2(x))  / 


1/30.  The  error 


(21) 


We  have  considered  30  realizations  of  <h  and  approximated  the  above  expec¬ 
tation  as  an  empirical  mean  over  these  30  realizations. 

In  Table  1,  we  compare  the  exact  solntion  of  (9)-(10)  (compnted 

using  a  hnite  element  method  with  a  hne  mesh  of  size  h  =  e/40  adapted  to 
the  small  scales  present  in  the  problem)  with  the  approximation 
obtained  using  the  approach  described  above  and  with  the  approximation 
'WmsFem(-,  1^)  obtained  using  the  standard  MsFEM  approach  (as  described 
in  Section  2.2),  with  the  recomputation  of  the  basis  functions  for  each  new 
realization  of  <h. 

Of  course,  the  computational  cost  to  obtain  many  realizations  of  wbll  is 
much  smaller  than  that  to  obtain  the  same  number  of  realizations  of  mmsFem- 
In  the  former  case,  and  in  contrast  to  the  latter  case,  we  do  not  have  to 
recompute  the  highly  oscillatory  basis  functions.  On  the  other  hand,  our 
approach  being  a  MsFEM  type  approach,  the  error  —  mmsFem  seems  to 
be  the  best  error  we  can  achieve.  It  is  thus  natural  to  compare  the  error  we 
obtain,  that  is  u'^  —  mbll,  with  that  “reference”  error. 


Example 

e(n*^,  mmsFem) 

e(n^MBLL) 

e(nMsFEM,  Wbll) 

Translation 

7.26% 

10.38% 

9.33% 

Rotation 

9.34% 

10.58% 

8.63% 

Table  1:  Errors  (21)  for  the  two  test  cases  considered. 


We  observe  that,  for  both  test  cases,  the  error  e(n'^,MBLL)  is  of  the  same 
order  as  the  reference  error  e(n'^,  mmsFem)-  Our  approach  thus  yields  an  ap¬ 
proximation  wbll  which  is  as  accurate  as  the  approximation  mmsFem  provided 
by  a  standard  MsFEM  approach,  for  a  much  smaller  computational  cost. 
Note  also  that  the  order  of  the  error  observed  here  (of  7  to  10  %  in  the 
norm)  is  the  standard  order  obtained  with  MsFEM  approaches  (see  e.g.  [2, 
Section  4.2,  Table  2]). 

These  hrst  numerical  results  are  encouraging.  Note  however  that  they 
both  correspond  to  the  specihc  case  where  V<F  is  piecewise  constant.  Dehnite 
conclusions  on  the  interest  of  the  approach  yet  need  to  be  obtained,  e.g.  using 
test  cases  with  more  complex  diffeomorphisms  <h. 
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4  Robustness  of  the  MsFEM  approach  to  a 
macroscopic  perturbation  of  the  diffusion 
coefficient 


[Work  expanded  in  [5].] 

This  section  is  devoted  to  a  preliminary  study  of  the  following  question. 
Assume  that  we  know  the  corrector  Wp  associated  to  a  periodic  coefficient 
^per  by  the  corrector  problem  (5).  Is  it  possible,  using  Wp,  to  approximate 
the  corrector  Wp  associated  to  a  macroscopic  perturbation  of  Aper?  Other¬ 
wise  formulated,  once  we  know  how  to  homogenize  (1)  with  the  coefficient 

A^{x)  =  Aper  is  it  possible  to  efficiently  homogenize  (1)  for  the  highly 

oscillatory  matrix  A^{x)  =  h{x)Ap^^  (~)’  +  ^per  Note 

that,  in  both  cases,  the  difference  between  A^{x)  and  A^{x)  only  comes  from 
a  function  b  that  has  no  small  scale  oscillation  (thus  the  terminology  “macro¬ 
scopic  perturbation”).  In  particular,  the  high  frequencies  present  in  A‘^{x) 
are  identical  to  the  high  frequencies  present  in  A^{x).  This  is  the  reason 
why  we  expect  that,  once  we  have  resolved  these  highly  oscillatory  modes  for 
A^(a;),  we  may  deduce  the  highly  oscillatory  modes  of  A^{x). 

Remark  5  The  same  question  may  be  asked  for  the  MsFEM  highly  oseil- 
latory  basis  funetions,  rather  than  the  periodie  eorreetor.  See  Seetion  f.2 
below. 

A  motivation  for  this  question  is  the  following  optimization  problem. 
Assume  that  we  model  our  material  with  some  coefficient  matrix  A^{x)  = 

(5)  +Mx).  that  we  kuow  that  the  highly  ose.llatory  eompouent  is  aeeu- 

rate,  and  that  we  want  to  optimize  on  the  macroscopic  component  b  in  order 
to  reproduce  with  this  model  some  known  results  (such  as  experimental  data) 
on  the  homogenized  behavior.  In  the  optimization  loop,  we  need  to  compute, 
for  each  new  trial  value  of  the  function  b,  the  homogenized  coefficient.  For 
the  sake  of  efficiency,  we  thus  need  to  perform  this  homogenization  procedure 
with  a  computational  cost  as  small  as  possible. 
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4.1  Multiplicative  pertubation 

As  a  first  step,  we  consider  a  multiplicative  perturbation.  The  coefficient 
A^{x)  reads 

A^{x)  =  b{x)Apev  ,  (22) 

where  &  is  a  scalar-valued  function.  The  problem  under  consideration  then 
reads 

— div  b{x)Ape^  j  Vu‘^{x)  =  f{x)  in  V,  u'^  =  0  on  dV.  (23) 

When  e  — 0,  the  function  converges  to  u*,  solution  to  the  homogenized 
problem 

—div  [A*(a;)VM*]  =  f{x)  in  V,  u*{x)  =  0  on  dV, 

where  the  homogenized  matrix  is  given  by  A*{x)  =  b{x)Ap^j.,  where  4*^,, 
is  defined  by  (4)-(5).  In  this  case,  the  corrector  associated  to  A‘^(x)  = 
b(x)Aper(x/e)  is  identical  to  the  corrector  associated  to  A^{x)  =  Aper{x/e). 

In  a  MsFEM  context,  we  have  the  same  type  of  result.  As  in  Section  2.2, 
introduce  the  PI  finite  element  space  Vh  =  Span(x°,  1  <i  <  L).  Compute 
the  highly  oscillatory  basis  function  xl  by  solving  (7).  These  functions  are 
therefore  independent  of  the  macroscopic  function  b.  We  next  introduce  the 
MsFEM  space  W/i  =  Span(xf,  1  <  f  <  T),  and  perform  a  Galerkin  approx¬ 
imation  of  (23)  using  the  space  Wh-  Because  of  the  specihc  structure  (22), 
such  an  approach  provides  an  approximation  of  solution  to  (23)  whose 
accuracy  is  essentially  independent  of  b.  For  each  new  trial  function  b,  we  do 
not  have  to  recompute  the  MsFEM  basis  functions. 

4.2  Additive  pertubation 

We  now  consider  an  additive  perturbation,  in  the  sense  that  the  coefficient 
A^{x)  reads 

A^[x)  =  b{x)  +  Aper  .  (24) 

We  do  not  assume  that  b  is  much  larger,  or  much  smaller,  than  Aper.  The 

ratio  —  is  of  the  order  of  one. 

II  Aper  II  L°° 
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For  the  sake  of  simplicity,  consider  first  the  case  when  6  is  a  constant 
function.  The  corrector  equation  associated  to  the  coefficient  (24)  then  reads 

f  -div[(Aper(|/)  +  b){p  +  Vwp{y))]  =  0  on 
\  Wp  is  Z'^-periodic. 

As  can  be  seen  on  the  one-dimensional  case,  there  is  no  relation  between  Wp 
and  the  corrector  Wp  associated  to  the  coefficient  A^(x)  =  Aper(x/e),  which 
solves  (5).  This  case  is  such  much  more  challenging  than  the  case  considered 
in  Section  4.1.  In  what  follows,  we  propose  a  numerical  approach,  based  on 
a  tensor-product  decomposition,  to  address  this  setting. 

Principle  In  the  sequel,  we  consider  coefficients  of  the  form 

A'^{x,  p)  =  Aq{x)  +  jj,b{x),  xgV,  /i  G  a,  (25) 

where  A  C  is  a  bounded  open  set.  The  coefficient  A^{x,  p)  is  thus  equal  to 
the  highly  oscillatory  coefficient  Aq{x),  up  to  the  addition  of  a  macroscopic, 
non  oscillatory  function  fj,b{x).  We  do  not  assume  that  Aq{x)  is  periodic,  and 
therefore  put  ourselves  in  the  MsFEM  context.  Our  aim  is  to  efficiently  com¬ 
pute  the  highly  oscillatory  basis  functions  xf(a;,  p),  solution,  in  each  element 
K,  to 


-div  [A^(a;,/i)Vxf(a;,p)]  =  0  in  K,  p)  =  X^i{x)  on  (9K,  (26) 

where  are  the  standard  PI  hnite  element  basis  functions.  In  turn,  the 
functions  xf(a;,  p)  will  be  used  to  perform  a  Galerkin  approximation  of  the 
problem 

—Aiv[A^{x^jj)Vu'^{x^p)]  =  f{x)  in  "D,  M^(-,p)  =  0  on  (9P,  (27) 

for  many  values  of  the  parameter  p  G  A. 

This  question  has  been  addressed  in  [15],  where  an  approach  based  on  the 
expansion  of  yf  (x,  p)  in  terms  of  a  Neumann  series  is  proposed.  Our  approach 
is  different,  and  consists  in  adapting  the  Proper  Generalized  Decomposition 
(PGD)  technique  [17,  12,  19,  18]  to  the  current  context.  More  precisely,  we 
are  going  to  approximate  yf,  function  of  the  two  variables  x  and  p,  as  a  sum 
of  products  of  a  function  depending  only  on  a;  by  a  function  depending  only 
on  fi. 
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We  now  proceed  in  details.  We  first  change  of  unknown  function  and 
define 

(28) 


£,K 

~  Xi  I K  X' 


i  Ik  Ik  ■ 


We  infer  from  (26)  that 


— div 


{Vx'i{x)  +  =0  in  K,  Vi’^{x,n)  =  0  on  dK. 

(29) 

The  advantage  of  considering  is  that  this  function  satishes  homogeneous 
boundary  conditions  on  dK.,  in  contrast  to  yf.  Our  approach  is  based  on  the 
assumption  that  nf’^(a;,/r)  writes  as  follows: 


£,K/ 


N 


(30) 


for  a  small  number  of  terms  N,  where  the  functions  fi  i— >■  Qijifi)  are  indepen¬ 
dent  of  X  and  the  functions  x  h- >  are  independent  of  /r.  Once  these 

functions  have  been  identihed,  computing  the  basis  functions  xf(-,/i)  for  any 
value  of  fi  just  amounts  to  evaluating  N  functions  of  fi  using  (28)  and  (30), 
rather  than  solving  the  partial  differential  equations  (26). 


Algorithm  The  functions  gij{fi)  and  fl'J^{x)  are  iteratively  dehned.  As¬ 
sume  that  they  have  been  built  for  any  j  <  k  —  1.  To  build  gi^k  and 
we  introduce  two  variational  formulations.  The  hrst  one  consists  in  hnding 
G  Hq(K)  solution  to 

WweH^,{K),  Ak{flf,w)  =  -Mw),  (31) 

with 

=  j^j^(^A%x,fi)Vf If  {x) -Vwix)^  glffi)dxdfi 

and 

J^kiw)  = 
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Note  that  it  is  natural  to  look  for  in  the  space  Hq(K),  since  satishes 
homogeneous  Dirichlet  conditions  on  dK.. 

The  second  variational  formulation  consists  in  hnding  ^  G  -h^(A)  such 
that 

Vh  e  L2(A),  h)  =  -7^fc(h),  (32) 

where 

^k{gi,k,  =  J  j  (^)  ■  9iAg)Kg)  dxd^i 

and 


TZkih)  — 

/k 


k-1 


i=i 


h(/i)  dx  djj,. 


Note  that  the  two  variational  formulations  (31)  and  (32)  are  coupled, 
as  they  both  involve  the  unknown  functions  and  gi^k-  However,  each 
unknown  function  only  depends  on  one  variable  {x  or  /i).  Solving  these  two 
coupled  problems  is  expected  (and  this  is  indeed  the  case)  to  be  easier  than 
solving  a  single  problem  on  a  function  depending  on  both  variables  x  and 
/i.  One  possibility  to  solve  (31)  and  (32)  is  to  use  the  following  iterative 
algorithm.  Let  g  be  the  accuracy  we  wish  to  reach  and  let  e  denote  the  error. 
We  proceed  as  follows: 

1.  Initialization:  set  e  =  1  and  gi^kig)  =  1/ so  that  ||5'i,fc||L2(A)  =  1. 

2.  Iterate: 


(a)  set  T  =  5fi,fc(/i); 

(b)  hnd  solution  to  (31); 

(c)  hnd  gi^k  solution  to  (32); 

(d)  multiply  the  function  gi^k  by  a  constant  such  that  its  norm  is 
equal  to  1:  gi^k  ^  gi,k/ \fj^k^ 


(e)  compute  the  difference  e 
old  iterate; 


[  {9^,k-Tr 
Ja 


between  the  new  and  the 
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(f)  if  e  >  T],  go  back  to  Step  2a. 

In  practice,  at  Steps  2b  and  2c,  both  problems  (31)  and  (32)  can  be  solved 
by  classical  methods.  For  instance,  we  can  discretize  the  bounded  domain 
A  C  and  use  a  hnite  element  method  to  solve  (32),  and  likewise  for  (31). 

Assume  now  that  the  functions  gi^k  and  have  been  computed  on  each 
element  K  and  for  each  k  <  N.  We  now  want  to  solve  (27),  for  some  value 
of  the  parameter  p  G  A.  We  introduce  the  MsFEM  space 

:=  Span{xf(a;,p)}. 

The  Galerlin  approximation  u%{x,iJ,)  G  of  the  solution  to  (27)  is 

dehned  as  the  solution  to 

vvewj;,  f  (VKfv(-,f,)v<(.,A)=  /  fv. 

Jv  Jv 


Numerical  illustration  We  work  in  dimension  two,  with  A^{xi,X2,  fi)  = 
a^{xi,X2i  jj)  Id,  where 

a^{xi,X2,  fi)  =  1  +  100  sin^  sin^  _l_  IQO  p,  {xi,X2)  G  M^, 

which  is  indeed  of  the  form  (25).  The  computational  domain  is  "D  =  (0, 1)^ 
and  the  parameter  domain  is  A  =  (0, 1).  We  set  e  =  0.05.  We  compute  the 
functions  Qi^k  and  as  explained  above,  and  evaluate  the  error 


e^(p) 


Vxf(-,/i)  - 

N 

Vx?  +  E 

7=1 

L2(K) 

1  (■)  a)  I  L^{K) 

(33) 


where  xf  (x,  p)  is  the  exact  basis  function  for  the  parameter  p,  which  solves  (26). 

We  have  worked  with  the  following  numerical  parameters.  The  open  set 
A  is  discretized  with  a  mesh  of  size  =  0.1.  To  solve  (26),  (31)  or  (32)  in 
practice,  each  hnite  element  K  is  discretized  with  a  mesh  of  size  h  =  77/30, 
where  H  =  diam(K).  We  take  H  =  1/7. 

In  Table  2,  we  show  the  error  (33)  as  a  function  of  N,  that  is  the  number 
of  terms  used  in  (30)  to  approximate  v/’^.  We  observe  a  fast  convergence 
of  the  error  with  respect  to  N.  Two  terms  in  (30)  are  actually  enough  to 
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N 

e^(/i  =  0.5) 

e^{p  =  0.67) 

1 

18.15% 

18.56% 

2 

3.78% 

4.18% 

3 

3.11% 

3.30% 

4 

2.82% 

3.13% 

Table  2:  Error  (33)  for  two  values  of  the  parameter  fi  (all  these  errors  cor¬ 
respond  to  the  same  choice  of  hnite  element  K;  similar  results  are  obtained 
for  a  different  choice). 

obtain  an  approximation  of  (and  therefore  of  xf)  with  an  error  smaller 
than  5  %. 

Obviously,  these  encouraging  results  are  only  preliminary.  More  tests  are 
needed  to  get  a  better  understanding  of  this  approach,  its  limitations  and 
the  regime  where  it  is  indeed  advantageous. 

5  Proposed  directions  of  research  for  an  ex¬ 
pected  renewed  funding 

If  BOARD  decides  to  renew  our  funding,  there  are  a  number  of  directions  of 
research  on  which  we  might  consider  proceeding  (subject  to  BOARD  approval 
of  course).  We  summarize  here  some  of  them. 

In  Section  4,  we  have  performed  a  preliminary  study  on  the  robustness 
of  the  MsBEM  approach  to  perturbations  that  are  non-oscillatory.  This 
question  is  actually  part  of  a  much  broader  question,  which  is  related  to 
inverse  problems  in  multiscale  science. 

A  first  remark  is  that  all  models  that  involve  a  random  parameter  require 
some  knowledge  on  the  distribution  of  this  random  parameter  (actually,  they 
most  often  require  a  complete  knowledge  of  that  distribution).  In  practice,  ac¬ 
cess  to  this  distribution  is  difhcult.  One  is  therefore  bound  to  assume  a  given 
form  (Gaussian,  . . . )  for  the  distribution  and  proceed  with  the  computation. 
A  question  of  major  practical  interest  is  to  a  posteriori  prove,  or  disprove 
the  validity  of  this  assumption.  Otherwise  stated,  tests  of  hypotheses  in  the 
context  of  engineering  problems  is  an  important  issue.  A  preliminary  step, 
before  trying  to  identify  the  distribution  of  the  random  parameters,  is  to  as¬ 
sume  a  specihc  form  of  this  distribution,  depending  on  a  few  quantities  (e.g. 
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assume  a  Gaussian  distribution  with  an  unknown  variance)  and  identify  these 
quantities.  To  solve  this  identihcation  problem,  it  is  important  to  be  able  to 
solve  efficiently  the  forward  problem  (given  the  microscopic  held  A,  compute 
the  macroscopic,  homogenized  behavior).  Efficient  methods  such  as  the  ones 
proposed  within  this  contract  are  then  of  paramount  importance. 

Another  remark  is  that  the  question  of  inverse  problems  in  materials  sci¬ 
ence  is  of  course  not  new.  However,  our  context  is  very  specihc,  owing  to 
the  fact  that  homogenization  acts  as  a  hlter.  Many  features  of  the  coefficient 

in  the  problem  (1)  (or  its  stochastic  version  (13))  are  hltered  out  by  the 
homogenization  procedure.  Several  helds  can  lead  to  the  same  homoge¬ 
nized  matrix  A*.  It  is  hence  hopeless  to  try  to  recover  the  held  A^  from  the 
sole  knowledge  of  A*,  or  of  properties  of  the  homogenized  material.  From 
A*,  one  can  only  expect  to  recover  the  class  of  microscopic  helds  that  cor¬ 
respond  to  this  homogenized  behavior.  This  class  probably  contains  many 
materials,  diherent  at  the  hne  scale,  but  equivalent  from  the  macroscopic 
standpoint.  All  of  these  are  thus  admissible,  if  the  only  information  we  have 
is  a  macroscopic  information. 
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Summary 

We  report  here  on  the  work  performed  during  the  second  year  (october  2010 
-  October  2011)  of  the  contract  FA  8655-10-C-4002  on  Multiscale  problems  in 
materials  science:  a  mathematical  approach  to  the  role  of  uncertainty. 

We  recall  that  the  bottom  line  of  our  work  is  to  develop  affordable  numer¬ 
ical  methods  in  the  context  of  stochastic  homogenization.  Many  partial  differ¬ 
ential  equations  of  materials  science  indeed  involve  highly  oscillatory  coeffi¬ 
cients  and  thus  small  length-scales.  Homogenization  theory  is  concerned  with 
the  derivation  of  averaged  equations  from  the  original  oscillatory  equations, 
and  their  treatment  by  adequate  numerical  approaches.  Stationary  ergodic 
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random  problems  (and  the  associated  stochastic  homogenization  theory)  are 
one  instance  for  modelling  uncertainty  of  continuous  media.  The  theoretical 
aspects  of  these  problems  are  now  well-understood,  in  general.  On  the  other 
hand,  the  numerical  aspects  have  received  less  attention  from  the  mathemat¬ 
ics  community.  Standard  methods  available  in  the  literature  often  lead  to 
very,  and  sometimes  prohibitively,  costly  computations. 

In  this  report,  we  first  focus  on  a  class  of  materials  of  moderate  difficulty 
but  of  signihcant  relevance,  that  of  random  materials  where  the  amount  of 
randomness  is  small.  They  can  be  considered  as  stochastic  perturbations  of 
deterministic  materials.  We  have  presented  in  the  previous  report  (see  [3]) 
a  possible  extension  of  the  well-known  Multiscale  Finite  Element  Method 
(MsFEM)  to  such  a  weakly  stochastic  setting,  along  with  detailed  numer¬ 
ical  results.  We  are  now  in  position  to  provide  a  complete  analysis  of  the 
approach,  extending  that  available  for  the  deterministic  setting. 

We  next  consider  a  different  weakly  stochastic  setting.  Rather  than  per¬ 
turbing  the  deterministic  material  by  frequent  but  small  random  amounts, 
we  consider  a  setting  in  which  the  deterministic  material  is  rarely  perturbed. 
However,  when  it  occurs,  the  perturbation  is  large.  Because  this  setting  is  a 
weakly  stochastic  setting,  the  workload  to  compute  the  homogenized  matrix 
is  already  smaller  than  in  generic  stochastic  homogenization.  We  show  here 
how  to  further  reduce  the  workload  by  using  a  Reduced  Basis  approach. 

We  hnally  turn  to  a  variant  of  stochastic  homogenization,  where  the  ran¬ 
domness  is  not  small,  and  describe  in  that  context  a  truncation  procedure  to 
compute,  in  practice,  an  approximation  of  the  homogenized  coefficient. 

The  works  described  below  have  been  performed  by  Claude  Le  Bris  (PI), 
Frederic  Legoll  (Co-PI)  and  Florian  Thomines  (second  year  Ph.D.  student). 

1  Introduction 

During  this  second  year  of  contract,  we  have  pursued  our  effort  on  developing 
affordable  numerical  methods  in  the  context  of  stochastic  homogenization. 

Many  partial  differential  equations  of  materials  science  indeed  involve 
highly  oscillatory  coefficients  and  small  length-scales.  Homogenization  the¬ 
ory  is  concerned  with  the  derivation  of  averaged  equations  from  the  original 
oscillatory  equations,  and  their  treatment  by  adequate  numerical  approaches. 
Stationary  ergodic  random  problems  are  one  of  the  most  famous  instances  of 
mathematical  uncertainty  of  continuous  media. 
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The  purpose  of  this  report  is  to  present  the  recent  progress  we  have  made 
on  this  topic,  with  the  aim  to  make  numerical  random  homogenization  more 
practical.  As  already  mentioned  in  the  previous  report,  because  we  cannot 
embrace  all  difficulties  at  once,  the  case  under  consideration  here  is  a  simple, 
linear,  scalar  second  order  elliptic  partial  differential  equation  in  divergence 
form,  for  which  a  sound  theoretical  groundwork  exists.  We  focus  here  on 
the  different  manners  the  problem  can  be  handled  from  the  computational 
viewpoint. 

We  begin,  in  Section  2,  with  a  brief  description  of  stochastic  homogeniza¬ 
tion,  the  only  purpose  of  which  is  the  consistency  of  this  report. 

As  pointed  ont  above,  random  homogenization  for  general  stochastic  ma¬ 
terials  is  very  costly.  Yet,  it  tnrns  ont  that  it  is  possible  to  identify  classes  of 
materials  of  moderate  difficnlty  bnt  of  signihcant  relevance,  where  stochastic 
homogenization  theory  and  practice  can  be  rednced  to  more  affordable,  less 
compntationally  demanding  problems.  These  materials  are  neither  periodic 
(because  such  an  oversimplifying  assnmption  is  rarely  met  in  practice),  nor 
fnlly  stochastic.  They  can  be  considered  as  an  intermediate  case,  that  of 
stochastic  perturbations  of  deterministic  (possibly  periodic)  materials.  Note 
that  many  practical  situations,  involving  actnal  materials  or  media,  can  be 
considered,  at  a  good  level  of  approximation,  as  perturbations  of  a  deter¬ 
ministic  (often  periodic)  setting  (see  e.g.  [15]).  In  this  report,  we  discuss 
two  different  weakly  stochastic  settings,  and  for  each  of  them,  we  present 
an  efficient  numerical  approach  to  handle  it.  First,  in  Section  3,  we  provide 
an  analysis  of  a  variant  of  the  Mnltiscale  Finite  Element  Method  (MsFEM), 
well  adapted  to  the  case  when  the  matrix  describing  the  properties  of  the 
material  is  the  snm  of  a  deterministic  term  and  a  small  random  term.  This 
variant  has  been  introdnced  in  the  previous  report  (see  [3]),  and  extensive 
numerical  tests  have  been  reported  there.  As  explained  below,  we  now  have 
a  complete  understanding  of  the  approach,  from  the  analysis  viewpoint.  We 
wish  to  emphasize  the  fact  that  considering  a  stochastic  pertnrbation  of  a 
deterministic  problem  and  handling  it  with  a  multiscale  technique  developed 
in  the  deterministic  setting  is  not  restricted  to  the  case  of  the  MsFEM.  Sim¬ 
ilar  entreprises  can  probably  be  nndertaken  in  other  settings,  snch  as  those 
proposed  in  [19]. 

In  Section  4,  we  tnrn  to  a  different  weakly  stochastic  setting,  introdnced 
by  the  PI  and  a  collaborator  of  his  in  [8,  9,  10].  This  model  is  well  snited 
for  representing  materials  with  rare,  bnt  non  small,  pertnrbations  with  re- 
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spect  to  a  deterministic  situation.  A  typical  example  is  a  composite  material 
embedding  fibers,  located,  say,  on  a  perfect,  periodic  lattice.  The  random 
perturbation  then  consists  in  deleting  some  fibers  (see  Figure  1  below).  This 
setting  is  a  weakly  stochastic  setting,  as  we  assume  that  such  an  accident 
occurs  very  rarely.  However,  it  is  clear  that  the  local  properties  (at  the  mi¬ 
croscopic  level)  of  the  material  are  significantly  changed  if  the  fiber  is  indeed 
deleted.  In  the  sequel,  we  show  that  the  Reduced  Basis  approach  can  be  used 
in  that  context  to  speed-up  the  computation  of  the  homogenized  coefficient. 

In  Section  5,  we  next  turn  to  a  non-weakly  stochastic  setting,  and  consider 
a  variant  of  stochastic  homogenization,  introduced  by  the  PI  and  co-workers 
in  [11,  12]  some  years  ago.  This  model  is  well  suited  to  represent  materials 
that  are  random  deformations  of  a  perfect  periodic  material.  A  typical  exam¬ 
ple  is,  again,  a  composite  material  with  fibers.  Fibers  are  all  identical,  they 
would  be  located  on  a  periodic  lattice  in  the  ideal  situation.  However  (for 
instance  as  a  consequence  of  the  manufacturing  process),  their  actual  posi¬ 
tions  are  now  random.  In  the  sequel,  we  present  and  analyze  a  procedure  to 
practically  approximate  the  homogenized  matrix. 

We  collect  in  Section  6  some  conclusions  about  the  work  performed  so 
far,  and  future  directions. 


2  Basics  of  stochastic  homogenization 

[Detailed  presentation  ean  he  read  in  [2,  3].] 


For  the  consistency  of  this  report  and  the  convenience  of  the  reader  not 
familiar  with  the  theory,  we  recall  here  some  groundwork  in  stochastic  ho¬ 
mogenization,  underlining  why  stochastic  homogenization  often  leads  to  ex¬ 
tremely  expensive  computations.  More  details  can  be  read  in  [3]  and  refer¬ 
ences  therein,  and  also  in  the  review  article  [2]  that  we  published. 

The  typical  random  homogenization  problem  writes 


— div 


A[—,uj]'Vu^\  =  f{x)  inV, 


u^{x)  =  0  on  dV, 


(1) 


where  A  is  a  bounded,  definite  positive,  stationary  (i.e.  statistically  homo¬ 
geneous)  random  matrix  (see  [3]).  In  this  framework,  it  is  well  known  that, 
as  e  — 0,  the  solution  u'^  to  (1)  converges  to  u*  solution  to 


— div[A*VM*]  =  f{x)  in  V,  u*{x)  =  0  on  dV,  (2) 
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where  the  homogenized  matrix  A*  reads 


l^%  =  E 


(ei  +  VWe,(y,-)f  A{y,-)  (ej +  Vw^j{y,-)) 


where  Q  =  (0,  and  where,  for  any  p  G  the  so-called  corrector  Wp  is 
the  (unique  up  to  the  addition  of  a  constant)  solution  to 


— div  [A  {y,uj)  {p +  'Vwp{y,uj))]  =  0  on 
Vwp  is  stationary,  E  (f  Vwp{y,  ■)  dy^  =  0. 


'Q 


(3) 


From  the  computational  viewpoint,  solving  (3)  is  challenging,  because  it  is 
posed  on  the  entire  space  R'^.  The  traditional  approach  is  to  truncate  (3) 
on  a  bounded  domain,  say  the  cube  Qn  =  {—N,Ny,  and  complement  it 
with  e.g.  periodic  boundary  conditions.  We  are  thus  left  with  solving  the 
truncated  corrector  problem 


J  — div  (74(-,c<;)  (p -I- Vtc^(-,a;)))  =  0  on  R'^, 

1  Wp{-,uj)  is  QAr-per iodic. 

In  turn,  the  homogenized  matrix  A*  is  approximated  by  the  matrix 


[■^*N]ij  (^) 


Although  A*  itself  is  a  deterministic  object,  its  practical  approximation 
A*^{u)  is  random.  It  is  only  in  the  limit  of  inhnitely  large  domains 
that  the  deterministic  value  is  attained.  Indeed,  as  shown  in  [17,  Theorem 
1],  we  have 

lim  =  A*.  (5) 

N—*oo 

Errors  between  A*j^{uj)  and  A*  are  due  to  (i)  the  truncation,  and  (ii)  the  fact 
that  the  truncated  problem  is  random  in  nature.  Because  of  the  truncation, 
E  [A^]  7^  A*.  At  hxed  N,  there  is  a  systematic  bias,  which  can  only  be 
reduced  by  taking  sufficiently  large  domains  Qv-  In  addition,  computing 
E  [A^]  is  also  expensive.  Indeed,  a  large  number  M  of  independent  realiza¬ 
tions  of  A^(a;)  should  be  considered  to  compute  an  empirical  mean,  in  the 
spirit  of  Monte  Carlo  methods.  It  is  only  in  the  limit  M  ^  oo  that  the  exact 
mean  E  [A^]  is  recovered. 

The  overall  computation  described  above,  that  involves  solving  several 
independent  realizations  of  (4)  on  presumably  large  a  domain  Qv,  is  thus 
very  expensive. 
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3  A  weakly-stochastic  MsFEM  approach 

[Work  expanded  in  [1,  4,  5].] 

Following  the  encouraging  numerical  results  reported  in  [3]  on  the  variant 
of  the  MsFEM  for  weakly  stochastic  settings,  we  have  pursued  our  efforts  and 
obtained  a  complete  analysis  of  the  proposed  approach,  that  we  describe  in 
the  sequel.  For  clarity,  we  begin  this  section  by  briefly  recalling  our  approach. 

3.1  The  proposed  approach 

We  consider  the  problem 

-div(A^(-,cc;)VM^(-,cc;))  =  /  vnV,  ^(•,cc;)  =  0  on  (9F>,  (6) 

where  G  is  a  random  matrix  satisfying  the  standard 

coercivity  and  boundedness  conditions.  In  contrast  to  (1),  we  do  not  assume 

that  A[^{x,u)  =  Ari  ,ci;^  for  a  fixed  stationary  matrix  A^j.  The  MsFEM 

approach  is  applicable  in  more  general  situations. 

We  suppose  that  A'[{x,u)  is  highly  oscillatory  in  both  its  deterministic 
and  stochastic  components,  and  that  it  is  a  perturbation  of  a  deterministic 
matrix,  in  the  sense  that 

A[^{x,uj)  =  Al{x) +riAl{x,uj),  (7) 

where  Aq  is  a  deterministic  matrix  and  r]  is  a  small  deterministic  parameter. 
This  model  may  be  well  suited  for  heterogeneous  materials  (or,  more  gener¬ 
ally,  media)  that,  although  not  periodic,  are  not  fully  stochastic,  in  the  sense 
that  they  may  be  considered  as  a  perturbation  of  a  deterministic  material. 

We  recall  that  the  MsFEM  approach  aims  at  approximating  the  solution 
of  (6)  by  performing  a  variational  approximation  of  the  problem  using  pre¬ 
computed  basis  functions  0f  that  are  adapted  to  the  problem.  The  main  idea 
of  our  proposed  approach  is  to  compute  a  set  of  deterministie  MsFEM  basis 
functions  0f  using  A^,  the  deterministic  part  of  A’^  in  the  expansion  (7),  and 
then  to  perform  Monte  Carlo  realizations  at  the  macroscale  level  using  a  set 
of  A4  realizations  of  the  random  matrix  ^ 

detailed  presentation).  Note  that,  for  each  of  these  realizations,  we  solve  the 
original  problem,  with  the  eomplete  matrix  and  not  only  its  deterministic 
part.  Only  the  basis  set  is  taken  deterministic. 
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The  deterministic  basis  functions  01  are  computed  only  once,  hence  the 
computational  saving  in  comparison  to  a  natural  adaptation  of  the  MsFEM 
to  the  stochastic  setting,  where,  for  each  realization  of  the  random  matrix 
new  basis  functions  are  computed  before  solving  the  macroscopic  problem. 

As  illustrated  by  the  numerical  tests  reported  in  [3,  4],  our  proposed 
approach  is  extremely  efficient  when  is  a  perturbation  of  Aq.  In  addition, 
the  small  parameter  r]  does  not  need  to  be  extremely  small  for  our  approach 
to  be  highly  competitive. 


3.2  Analysis 

We  now  turn  to  the  analysis  of  our  approach.  We  recall  that,  in  the  deter¬ 
ministic  setting,  a  classical  context  for  proving  convergence  of  the  MsFEM 
approach  (see  [21])  is  the  case  when,  in  the  reference  highly  oscillatory  prob¬ 
lem 

-div(AVu")  =  /  in  T>,  u"  =  0  on  dV,  (8) 

the  matrix  reads  A^{x)  =  Ap^r  ^  for  a  fixed  periodic  matrix  Ap^r-  Likewise, 
to  be  able  to  perform  our  theoretical  analysis  in  the  stochastic  setting,  we 
assume  that  Ap{x,uj)  =  Ap  for  a  hxed  stationary  random  matrix  Ap, 

although,  we  repeat  it,  the  approach  can  be  used  in  practice  for  more  general 
cases.  The  problem  (6)  then  admits  a  homogenized  limit  when  e  vanishes. 

Our  proof  follows  the  same  lines  as  that  in  the  deterministic  setting,  which 
we  now  briefly  review.  The  MsEEM  is  a  Galerkin  approximation,  the  error 
of  which  is  then  estimated  using  the  Cea  lemma: 


-  mmIIhu©)  <  o  inf  Wu^-vhWh^v),  (9) 

Vh&Wh 

where  u'^  is  the  solution  to  the  reference  deterministic  highly  oscillatory  prob¬ 
lem  (8),  um  is  the  MsFEM  solution,  Wh  =  Span((;/)f)  is  the  MsFEM  basis  set, 
and  C*  is  a  constant  independent  of  the  small  length-scale  e  present  in  A^ 
and  of  the  macroscopic  mesh-size  h.  Taking  advantage  of  the  homogenization 
setting,  we  introduce  the  two-scale  expansion 


=  u* 


w 


ei 


i=l 


du* 

dxi 
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of  where  u*  is  the  homogenized  solution  and  w^.  is  the  periodic  corrector 
associated  to  G  We  deduce  from  (9)  that 

\\u^ -umWhHv)  <C  i\\u^ -v^\\Him)+  inf  \\v^  -  VklluHv) 

\  vh&Wh 

The  hrst  term  in  the  above  right-hand  side  is  estimated  using  standard  ho¬ 
mogenization  results  on  the  rate  of  convergence  of  —  u^.  To  estimate  the 
second  term,  one  considers  a  suitably  chosen  element  Vh  G  W/i,  for  which 
~  Vh\\H^  can  be  estimated  directly.  The  main  idea  is  that  the  highly 
oscillating  part  of  can  be  well  approached  by  an  element  in  W/j,  since, 
by  construction,  the  highly  oscillatory  basis  functions  01  are  dehned  by  a 
problem  similar  to  the  corrector  problem,  and  thus  encode  the  same  highly 
oscillatory  behavior  as  that  present  in  the  correctors  w^..  We  are  thus  left 
with  approximating  the  slowly  varying  components  of  for  which  standard 
FEM  estimates  are  used. 

Following  the  same  strategy  in  our  stochastic  setting,  we  estimate  the 
distance  between  the  solution  to  the  reference  stochastic  problem  (6)-(7) 
and  the  weakly  stochastic  MsFEM  solution  us  as 

vh&yyh  / 

We  observe  that  a  key  ingredient  for  the  proof  is  the  rate  of  convergence  of 
the  difference  between  the  reference  solution  and  its  two-scale  expansion 
Such  a  result  is  classical  in  periodic  homogenization,  but,  to  the  best  of 
our  knowledge,  open  in  the  general  stationary  case  (in  dimensions  higher  than 
one).  One  only  knows  that  — vanishes  (in  some  appropriate  norm)  when 
e  — 0.  However,  in  the  particular  case  when  is  only  weakly  stochastic,  we 
have  shown  in  [5,  Theorem  2]  such  a  result,  useful  to  control  the  hrst  term 
in  (10): 


<  C'  (v^  +  r/x/£ln(l/e)  +  , 

where  C*  is  a  constant  independent  of  e  and  rj.  This  result  relies  on  asymptotic 
properties  of  the  Green  function  of  the  operator  L  =  — div  [HperV-],  a  topic 
of  independent  interest  which  has  been  investigated  in  [1]. 
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Hence,  exploiting  the  specificity  of  our  weakly  stochastic  setting,  we  have 
estimated  the  error  given  by  our  approach  as  (see  our  main  result.  Theorem  10 
in  [4]): 


-  Us 


<  c 


(  £  t  £  \ 

(  A/i+ h  + - +  ?7  \\i{N {h))  +  rj  +  rf‘C{r]) 


•> 


where  C  is  a  constant  independent  of  e,  h  and  77,  C  is  a  bounded  function  as 
f]  goes  to  0,  N{h)  is  the  number  of  elements  in  the  mesh  (roughly  of  order 
h-d  dimension  d),  and  ||  ■  is  a  broken  norm,  defined  by 


u 


Hi 


1/2 


where,  in  the  above  sum,  K  is  any  element  of  the  coarse  mesh  T^. 


Remark  1  As  is  often  the  ease  in  the  deterministic  MsFEM,  we  use  in  [4] 
the  oversampling  technique,  which  is  known  to  improve  the  accuracy  of  the 
numerical  results.  Consequently,  the  basis  functions  do  not  belong  to 
hence  the  use  of  a  broken  norm  in  the  above  estimate.  We  refer 
to  [4]  for  more  details. 

It  is  worth  noticing  that,  when  17  =  0  in  (7),  our  approach  reduces  to  the 
standard  deterministic  MsFEM  (with  oversampling),  and  the  above  estimate 
then  agrees  with  those  proved  in  [21], 


4  Reduced  Basis  approach  in  a  weakly  stochas¬ 
tic  homogenization  setting 

[Work  expanded  in  [6].] 

4.1  Summary  of  previous  works 

In  the  previous  works  [8,  9,  10],  not  funded  by  BOARD,  the  PI  and  a  collab¬ 
orator  of  his  introduced  the  following  weakly  stochastic  case. 

Consider  the  highly  oscillatory  problem  (1),  where  the  matrix  A  reads 

A{x,U}'^  Aper(27)  -|-  6,^ (x,  Co’)C*per (s^)  (H) 
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where  Aper  and  Cper  are  two  periodic  matrices,  and 

brj{x,u)  =  lQ+k{x)B^{u) 

where  are  i.i.d.  scalar  random  variables,  sharing  the  following  law: 

=  1  with  probability  77,  and  =  0  with  probability  1  —  77.  In  the  sequel, 
77  is  a  small  parameter,  so  that  A  =  Ap^,.  “most  of  the  time” .  We  hence  see 
that  the  perturbation  introduced  by  brj{x,uj)Cpe^{x)  in  (11)  is  rare.  On  the 
other  hand,  since  Aper  +  Cper  is  very  different  from  Aper,  the  perturbation, 
when  it  occurs,  is  large.  See  Figure  1  for  some  illustration. 


Figure  1:  From  left  to  right;  perfect  material  (modelled  with  Ape,-),  material 
with  one  defect  and  two  defects. 


As  explained  in  Section  2,  we  approximate  A*  using  the  standard  trunca¬ 
tion  method  for  the  corrector  problem  (see  (4)).  By  enumerating  all  possible 
realizations  of  A{x,  u)  on  Qn,  we  obtain  an  expansion  of  E  [A)(r]  in  powers  of 
77  (see  [9,  10,  2]): 


where 


A^e,  = 


Al’^ei 


E  [41*^]  =  a;,,  +  77A*’^  +  r^^AY  +  •  •  •  , 


(12) 


/  A,{Vwlf  +  e,)  -  /  Aper(V<  +  e,), 

JQn  •JQn 

\Y{I  A2’"(V7n2f’^-fei) -2  f  AiiVwlf  +  ei) 

g— 1  x-'Qn  •'Qn 


+  /  Aper(V7ng.  -f  Cj)  )  , 

J  Qn 


(13) 
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where  is  the  corrector  associated  to  (perfect  material),  solution  to 

-div  [Aper  {p  +  VtCp)]  =  0,  Wp  is  Q-periodic,  (14) 

and  is  the  corrector  associated  to  Ai  =  A^^^  +  IqC^er  (material  with 

one  defect): 

—div  [Ai  {p  +  VWp’^)]  =  0,  Wp’^  is  (^Ar-periodic.  (15) 

In  turn,  is  the  corrector  associated  to  AI’^  =  Aper  +  IgCper  +  Ig+sC'per 

(material  with  two  defects,  located  in  Q  and  Q  +  s) : 

—div  [^2’*  [p  +  VtCp  =  0,  Wp’^’^  is  QAr-periodic.  (16) 

Note  that  the  periodic  boundary  conditions  in  (15)-(16)  allow  us  to  assume, 
without  loss  of  generality,  that  the  hrst  defect  is  located  in  the  cell  Q. 

It  has  been  shown  numerically  in  [9,  10]  that  the  expansion  (12)  is  not 
only  valid  in  the  asymptotic  regime  p  1,  but  also  for  practical  small  values 
of  p.  In  some  cases,  the  expansion  is  even  valid  for  values  of  p  as  large  as 
0.5,  in  which  case  the  random  variables  take  value  0  and  1  with  eqnal 
probability. 

Note  that  the  compntation  of  ^2’'^,  when  necessary,  reqnires  to  solve  the 
corrector  problems  (16)  for  any  value  of  s  (the  position  of  the  second  defect). 
In  the  seqnel,  we  propose  to  use  a  Reduced  Basis  approach  to  solve  these 
Nd  _  I  problems,  that  are  parameterized  by  s. 

All  the  results  of  this  section  are  illustrated  with  the  same  two-dimensional 
numerical  example,  that  we  now  introduce.  We  take 

Aper(a:,  y)  =  20  Id2  -h  100  ^  Iq+kix,  y)  sin^(7ra;)  sin^(7r|/)  Id2 

k&I? 

and 

Cpe,{x,y)  =  -100  ^  lQ+k{x,y)  sin^(7ra;)  sin^(7r|/)  Id2. 

k&I? 

In  line  with  Figure  1,  this  test  case  represents  a  material  with  constant 
properties,  reinforced  by  a  periodic  lattice  of  circular  inclusions.  Loosely 
speaking,  the  pertnrbation  consists  in  randomly  eliminating  some  hbers.  See 
Figure  2  for  a  particular  realization  of  the  material. 

In  the  seqnel,  we  focus  on  the  hrst  entry  [A)(r]ii  homogenized  ma¬ 

trix.  We  thus  set  p  =  Ci  in  (14),  (15)  and  (16).  These  corrector  problems 
are  numerically  solved  using  a  mesh  of  size  h  =  1/10.  Qualitatively  similar 
conclusions  are  obtained  with  the  other  entries. 
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Figure  2:  Left:  the  perfect  (periodic)  material.  Right:  a  realization  of  the 
material  with  some  defects. 

4.2  Difficulties 

In  [6],  we  propose  to  use  a  Reduced  Basis  approach  to  speed-up  the  compu¬ 
tation  of  the  family  of  problems  (16)  parameterized  by  s,  the  weak  form  of 
which  is 

u;  s)  =  bp{v),  (17) 

where 


The  Reduced  Basis  approach  (see  [13]  for  a  presentation  of  the  method  in 
the  stochastic  case)  can  be  understood  as  a  way  to  approximate  the  set  of 
functions 

g  ■=  <  s<N‘^-l} 

by  an  element  in  the  space 

Xm  =  SpanlWp’’^’"’'^,  1  <  m  <  M},  (18) 

for  some  well-chosen  values  of  Sm,  1  <  — 1.  This  approach  is  efficient 

if  we  can  choose  a  small  value  for  the  dimension  M  of  Xm,  while  maintaining 
accuracy. 

Once  Xm  has  been  defined,  we  approximate  the  solution  of  (17)  us¬ 
ing  a  standard  Galerkin  approximation  on  Xm-  we  approximate  by 

^2,s,N,M  g  solution  to 

Wvm  e  Xm,  a(wJ^’^’^,nM;s)  =  bp{vM)-  (19) 
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The  construction  of  is  performed  beforehand  by  computing  the  solution 
of  (17)  for  some  values  of  the  parameter  s.  The  set  of  appropriate  param¬ 
eters  {sm} i<m<M  selected  using  a  standard  procedure,  called  the  Greedy 
procedure.  The  Reduced  Basis  approach  is  expected  to  be  efficient  because 
M  is  small  (the  problem  (19)  is  hence  posed  in  a  low- dimensional  space,  and 
thus  easy  to  solve)  and  because  there  are  —  1  problems  (17)  to  be  solved. 

Blindly  applied  to  the  family  of  problems  (17),  the  Reduced  Basis  ap¬ 
proach  is  not  efficient,  because  we  would  have  to  take  M  of  the  order  of 
iV*^.  The  reason  is  the  following.  We  first  compute  the  solution  of  the  prob¬ 
lem  (17)  for  all  the  different  values  of  the  parameter  s,  and  then  proceed 
with  a  Proper  Orthogonal  Decomposition  of  the  family  for 

the  scalar  product.  We  then  look  at  the  decay  of  the  spectrum,  which  is 
shown  on  Figure  3.  We  observe  no  decay  of  the  spectrum,  which  indicates 
that  the  functions  are  linearly  independent.  There  is  no  good  struc¬ 

ture  in  that  family,  and  therefore  we  would  need  to  choose  M  of  the  order 
of  to  accurately  approximate  the  functions  in  £  by  an  element  in  Xm- 
There  is  thus  no  speed-up.  The  appropriate  way  of  applying  the  Reduced 
Basis  approach  is  described  in  Section  4.3.1. 


Figure  3:  POD  (using  the  scalar  product)  of  the  family  (^p 
for  N  =  11. 

A  second  difficulty  is  that  the  Reduced  Basis  approach  requires  an  a  pos¬ 
teriori  error  estimator,  both  to  select  the  appropriate  values  of  Sm  in  (18), 
and  to  assess  the  quality  of  the  final  outcome.  In  our  particular  context, 
because  of  the  specificities  of  our  problem,  the  computation  of  the  classi¬ 
cal  error  estimator  turns  out  to  be  prohibitively  expensive.  We  describe  in 
Section  4.3.2  a  way  to  circumvent  this  difficulty. 
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4.3  Adjustment  of  the  Reduced  Basis  approach 

In  this  section,  we  show  how  to  address  the  two  difficulties  underlined  above. 


4.3.1  Building  a  family  of  functions  with  a  good  structure 

We  have  seen  above  that  there  is  no  good  structure  in  the  family  of  functions 
.  In  this  section,  we  rewrite  these  functions  as  linear  combinations 
of  the  functions  and  dehned  below.  The 

interest  is  that  these  two  FaFter  families  have  a  good  structure. 

To  build  these  good  families,  we  begin  by  introducing 


P 

^2,s,N 


=  - 
P 

= 


w. 


pi 


w, 


1,N 


W, 


1,N( 


s)  +  wj. 


where  Wp,  Wp’^  and  Wp’"^’^  are  dehned  by  (14),  (15)  and  (16),  respectively. 
Heuristically,  this  amounts  to  subtracting  the  appropriate  reference  function 
from  the  correctors  Wp’^  and  Wp’^’^ .  One  can  show  that  Wp'^  and  go 

to  0  at  inhnity,  and  are  hence  essentially  supported  in  a  compact  domain,  in 
contrast  to  Wp’^  and  Wp'"^’^.  We  see  that  Wp'^  is  solution  to 


-  div  (AiVwJ’^)  =  div(lQOper(VwJ  +p)), 
Qn  —  periodic. 


and  Wp'^’^  is  solution  to 

f  -  div  =  div  (liQ+^C'perVwJ’^)  +  div  {iqCperVwl’^ {■  -  s))  , 

1  Wp^’^  Qn  —  periodic. 


We  next  use  the  linearity  of  (20),  and  rewrite  as 


where  Wp'^’^  solves 

I  -  div  =  div  ((liQ+^OperVluJ’^)  , 

Qn  —  periodic. 
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and  solves 

r  -  div  ’^)  =  div  (iQCperVw;’^(-  -  5)) , 

\wp'^’^  Qv  —  periodic. 

Due  to  the  periodic  boundary  conditions,  we  have  Wp’^’^  =  Wp~'^’^{-  —  s). 
After  tedious  but  straightforward  computations,  and  assuming  for  the  sake  of 
simplicity  that  Aper  is  symmetric,  we  recast  the  matrix  defined  by  (13) 
as 


N'^-l  „ 

■4A)  =  E/  (le  +  lo+.)(V'<+ejfCp.,V5iy'' 

Jqn 

+  [  {Vwl+e,fC,,,Vwlf{--s). 
Jq 

Only  the  functions  Wp’^  and  needed  to  evaluate 

A  similar  formula  holds  if  Aper  is  not  symmetric,  up  to  slight  modihcations. 

We  now  numerically  test  for  the  linear  independence,  or  dependence,  of 

with  iV  =  11.  There 


the  functions  [wp’^’^) 


l<s<N^-l 


and  (wp’^’^) 


l<s<N^-V 


are  N^  —  1  =  120  functions  is  each  family.  We  compute  the  POD  of  the  family 

of  the  family  plot  in  Figure  4  the 

120  POD  eigenvalues  in  decreasing  order. 

In  contrast  to  Figure  3,  we  now  observe  a  decay  in  the  eigenvalues.  There 


is  a  good  structure  in  both  families 


1<S<V2-1 


and 


i<s<Af2-r 


Therefore,  much  fewer  vectors  are  needed  to  approximate  these  families  than 
for  approximating  the  family  considered  in  Section  4.2. 

Otherwise  stated,  we  expect  that  the  dimension  of  the  Reduced  Basis  ap¬ 
proximation  space  that  we  will  need  is  low,  in  contrast  to  the  situation  of 
Section  4.2.  Note  also  that  the  decay  is  somewhat  stronger  for  the  family 


W’'’'^)i<s<v2-i  ^han  for 


l<s<Af2-l' 


4.3.2  Modifying  the  error  estimator 

We  now  address  the  second  difficulty  identihed  at  the  end  of  Section  4.2, 
related  to  the  classical  a  posteriori  estimator,  and  suggest  an  alternative 
estimator,  better  suited  to  our  context  here. 
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Figure  4:  POD  (using  the  scalar  product)  of  the  family 
(left)  and  of  the  family  (right),  for  N  =  11. 


The  idea  is  to  use  the  norm  of  the  discrete  residual,  rather  than  the  clas¬ 
sical  estimator,  based  on  the  norm  of  the  dual  of  the  residual.  Let  us  denote 
Wh  '■=  span((;/)j)  the  finite  element  space  used  to  solve  the  problem  (16),  in 
the  reference  approach.  We  dehne  the  ith  component  of  the  discrete  residual 
as 


GT{s)  :=a(n;f’^’™,0p5)-6,(0.), 

where  jg  approximation  of  the  solution  to  (16)  computed  using 

the  approximation  space  (see  (19)).  In  the  sequel,  we  work  with  the  error 
estimator  _ 


A^(s):=  ||G”^(s)||,2 


The  classical  error  estimator  is  denoted  Am(<s). 


4.4  Numerical  results 

We  collect  here  numerical  results  that  demonstrate  the  efficiency  of  our  ap¬ 
proach. 

We  hrst  compute  the  relative  error  on  between  the  reference  value  (13) 
and  its  approximation  obtained  using  the  Reduced  Basis  approach,  for  differ¬ 
ent  values  of  the  dimension  M  of  the  approximation  space  Xm-  In  Figure  5, 
we  plot  this  indicator  for  two  components  of  the  matrix  .  We  observe 
that,  for  both  estimators  Am{s)  and  Am{s),  the  accuracy  improves  rapidly 
with  M.  Choosing  a  value  for  M  of  the  order  of  20  is  enough  to  obtain  a  good 
accuracy.  Since  this  value  is  much  smaller  that  the  total  number  of  functions 
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to  compute,  which  is  here  iV^  —  1  =  120,  we  obtain  a  good  computational 
speed-up. 


Figure  5:  Relative  error  on  as  a  function  of  the  dimension  M  of  the 
Reduced  Basis  space  Xm-  Left:  Right:  {A2^)i2. 


Next,  in  Figure  6,  we  plot  the  indicator 


e(m) 


max 

1<S<V2-1 


||Vw;2.s,iV 


Y7„,,2,s,V,m||  „ 
VWp’  ||l2(q^) 


I  ^-7  2,s,iV  I 

I  ^7  'Wp  I 


L‘^{Qn) 


(21) 


for  the  two  Reduced  Basis  approaches  (using  the  error  estimators  or 
Am).  In  both  cases,  we  observe  a  signihcant  decay  of  e(m)  as  a  function  of 
m.  Starting  from  an  error  of  100%  with  m  =  1  (the  Reduced  Basis  space  is  of 
dimension  one),  we  see  that  the  error  decreases  down  to  2%  for  M  20.  As 
above,  we  hence  see  that  a  small  value  of  M  is  sufficient  to  ensure  accuracy, 
hence  the  computational  speed-up. 

On  Figure  7,  we  show  in  which  order  the  relevant  cells  (i.e.  parameters 
Sm)  are  selected  by  the  Greedy  procedure.  We  compare  in  this  matter  the 
two  error  estimators.  We  observe  that,  although  the  actual  order  may  vary, 
the  two  estimators  provide  the  same  set  of  parameters  both  for 

the  eight  and  the  twenty  best  parameters. 


5  A  variant  of  stochastic  homogenization 

[Work  expanded  in  [7].] 

As  announced  in  the  previous  report  [3],  we  have  considered  a  variant 
of  the  classical  setting  of  stochastic  homogenization,  originally  introduced  a 
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Figure  6:  Relative  error  (21)  on  the  two-defects  correctors  .  The  mag¬ 
nitude  of  the  residual  error  for  large  M  (of  the  order  of  2  %)  is  close  to  the 
Finite  Element  error. 
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Figure  7:  Ordering  of  the  cells  selected  by  the  Greedy  procedure,  for  the  two 
error  estimators,  Am{s)  (left)  and  Am(s)  (right). 


few  years  ago  in  [11,  12].  The  equation  under  consideration  is  (1),  where 
the  matrix  A  is  the  composition  of  a  periodic  matrix  Aper  with  a  stochastic 
diffeomorphism  <F; 


A 


:=  A 


per 


(22) 


We  assume  that,  almost  surely,  the  map  <h(-,a;)  is  a  well-behaved  diffeo¬ 
morphism  (in  the  sense  that  Esslnf^^f^  (det(V<F(x,  a;)))  =  ly  >  0  and 
EssSup^gj^ 3,gRd  |V<F(x, a;)|  =  M  <  -|-oo),  and  that  it  satisfies 

V<h  is  stationary.  (23) 


Formally,  such  a  setting  is  well  suited  to  model  materials  that  are  periodic, 
in  a  given  reference  configuration.  The  latter  is  only  known  up  to  a  cer- 
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tain  randomness.  Materials  we  have  in  mind  are  ideally  periodic  materials, 
where  some  random  deformation  has  been  introduced,  for  instance  during 
the  manufacturing  process.  Assumption  (23)  means  that  V<h  is  statistically 
homogeneous,  i.e.  the  randomness  is  the  same  anywhere  in  the  material. 

The  problem  (l)-(22)  admits  a  homogenized  limit  when  e  vanishes.  It 
is  indeed  shown  in  [11]  that,  under  the  above  assumptions,  solu¬ 

tion  to  (l)-(22)  converges  as  £  goes  to  0  to  m*,  solution  to  the  homogenized 
problem  (2).  The  homogenized  matrix  A*  is  given  by,  for  any  1  <  i,  j  <  d, 


A*j  =  det  (  E 


'Q 

E 


-1 


X 


^^iy,-)dy 

/ 

ef  Aper  (y,  ■))  (ej  +  Vwe,  (l/,  ■))  dy  )  ,  (24) 


where  Q  =  (0, 1)*^  and  where,  for  any  p  G  W^,  Wp  solves  the  corrector  problem 


-div  [Aper  ($  ^{y,uj))  {p  +  Vwp{y,u))]  =  0  in 
Wp{y,u)  =  Wp{^~^{y,u),uj),  Vwp  is  stationary. 


E 


Vwp{y,  ■)dy  =  0. 


(25) 


The  hrst  question  we  address  in  [7]  is  to  precisely  understand  the  behavior 
of  u^{x,uj)  —  u*{x)  when  e  goes  to  0,  where  is  the  solution  to  the  highly 
oscillating  problem  and  u*  the  solution  to  the  homogenized  problem.  In  the 

1-11  T  •  1  r  u'^(x,uj)  —  u*(x) 

classical  random  ergodic  setting,  the  convergence  of  - -= -  to  a 

Gaussian  process  has  been  shown  in  [16,  18]  for  the  one-dimensional  case. 
We  show  in  [7]  a  similar  result,  for  the  variant  we  consider  here. 

We  have  also  investigated  the  question  of  how  to  approximate  the  homog¬ 
enized  matrix  (24)  in  practice.  Note  indeed  that  the  corrector  problem  (25) 
is  posed  on  the  entire  space  and  thus  challenging  to  solve.  In  the  classical 
context,  described  in  Section  2,  the  corrector  problem  (3)  is  also  posed  on 
and  a  standard  procedure  is  to  solve  the  corrector  problem  on  a  truncated 
domain  Qv  =  {—N,NY  (see  (4)).  The  convergence  of  the  procedure  when 
N  ^  oo  is  given  by  [17,  Theorem  1].  In  [7],  we  perform  a  similar  analysis  in 
the  context  of  (24)- (25).  In  the  following,  we  hrst  describe  a  numerical  strat¬ 
egy  to  approximate  the  homogenized  matrix  A*,  which  was  hrst  introduced 
in  [14].  We  next  provide  an  example  of  numerical  results,  before  turning  to 
the  analysis  of  the  approach. 
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Derivation  of  the  truncated  problem  The  weak  formulation  of  the 
corrector  problem  (25)  reads  as  follows  (see  [11]):  for  all  ip  stationary,  we 
have 

¥.(  j  iyipiy))'^ Aper{^~^{y,uj)){j)  +  Vwp{y,uj))d^  =0, 


where  ip  =  ip  o  ^  The  above  expression  can  be  rewritten,  after  a  change 
of  variables,  as 


E 


det(V<h)  (V$)-'^klper  (p+  (V<h)-Vwp) 


0. 


As  Ip,  V<h,  Aper  and  Vwp  are  stationary,  this  equivalently  reads,  because  of 
the  ergodic  theorem, 

lim  ^  [  det(V$)  (vipY  (V$)-%er  {p  +  {V^y^Vwp)  =  0. 

N^oo  \Qj^\  \  J 

At  hxed  TV,  we  now  dehne  the  approximate  corrector  as  the  Qv-periodic 
function  satisfying 


det(V<h) 


(v^)^  (V0)-%er  {p  +  (V0)-VO 


0 


for  all  Ip  Qv-periodic,  or,  equivalently, 

r  -div  [det(V$)(V0)-'^Aper  {p  +  (V0)-iV^if )]  =  0, 
1  Wp{-,uj)  is  Qv-periodic. 


(26) 


In  turn,  we  approximate  the  homogenized  matrix  A*  dehned  by  (24)  by,  for 
any  1  <  qj  <  d. 


[^v(^)]p'  =  det  V0(-,o;)^  x 

T^l  [  e^jAper{(p~\y:Uj))  (ej  +  Vw^{y,uj))  dy  (27) 

where  Wp{y,uj)  =  Wp  {(p~^ {y ,  u) ,  u) .  Note  that,  as  is  standard  when  ap¬ 
proximating  the  homogenized  matrix  of  stochastic  elliptic  problems,  the  ap¬ 
proximation  A*j,f{u)  is  a  random  matrix,  even  though  the  exact  homogenized 
matrix  (24)  is  deterministic.  This  is  a  by-product  of  working  on  the  truncated 
domain  Qv  rather  than  W^. 
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Numerical  illustration  Some  numerical  tests  following  our  strategy  (26)- 
(27)  are  reported  in  [2,  Section  3.2],  We  reproduce  here  some  of  them. 

Consider  the  following  two  dimensional  test-case:  we  give  ourselves  two 
families  {Xk)k£z  and  (yfc)fcez  of  i-i.d.  random  variables,  all  sharing  the  uni¬ 
form  law  U{[a,b]),  with  a  =  —2.25  and  b  =  5.75.  We  then  consider  the 
diffeomorphism  <h(a:,a;)  =  Qx  +  with  x  =  {x\,X2)  G  ^{x,uj)  = 

{i)x{xi,uj),iiY{x2,u])),  where 

f  [fc,A:-|-l[(^l) 

fcez 

and  likewise  for  ipy  The  periodic  matrix  Aper  is  defined,  for  all  x  E  Q,hy 

Aper{x)  =  aper{x)  Id2,  aper(xi,  X2)  =  (3  +  {a  -  (3)  sin^(7ra;i)  sin^(7rx2), 

with  a  =  10  and  (3=1.  The  results  obtained  with  (26)-(27)  are  shown  on 
Fig.  8.  We  indeed  observe  the  convergence  of  our  approximation  as  N  ^  oo. 


i)x{xYUj)  =  ^ 


Figure  8:  For  any  iV,  following  (26)-(27),  we  compute  several  realizations  of 
A^(a;),  and  build  from  these  a  confidence  interval  for  E  [(^^)]^i]. 


Convergence  of  the  numerical  strategy  In  [7],  we  study  the  conver¬ 
gence  of  our  approach,  as  N  ^  oo.  We  obtained  the  following  result,  which 
generalizes  [17,  Theorem  1]  (see  (5)  above)  to  the  variant  we  consider  here. 
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Theorem  1  Let  ^  be  a  dijfeomorphism  that  satisfies  (23),  and  Aper  be  a 
periodic,  coercive  and  bounded  matrix.  Then  the  random  matrix  de¬ 

fined  by  (27)  converges  almost  surely  when  N  ^  oo  to  the  deterministic 
homogenized  matrix  A*  defined  by  (24). 

6  Conclusions  and  agenda  for  the  third  year 
of  contract 

We  summarize  here  the  directions  of  research  we  wish  to  pursue  during  the 
third  year  of  contract. 

A  hrst  track  is  to  develop  a  MsFEM-type  approach  for  random  materials 
modelled  by  (l)-(22)  (see  Section  5).  Our  idea  relies  on  approximating  the 
solution  to  the  corrector  problem  (25)  by  the  periodic  corrector  (associated 
to  the  periodic  matrix  Aper)  composed  with  the  random  diffeomorphism  <h. 
The  advantage  of  this  approximation  is  that,  for  any  new  realization  of  <h,  we 
do  not  have  to  recompute  the  corrector.  We  have  developed  a  MsFEM-type 
approach  using  this  approximation,  and  have  performed  some  preliminary 
numerical  experiments,  that  yield  encouraging  results.  More  comprehensive 
tests  are  yet  needed  to  better  understand  the  capabilities  of  this  approach. 

A  second  direction  concerns  the  use  of  Reduced  Basis  methods  in  a  multi¬ 
scale  context.  We  have  seen  in  Section  4  that  the  Reduced  Basis  approach  can 
be  used  to  speed-up  the  computation  of  the  corrector  problems.  In  the  next 
future,  we  would  like  to  apply  this  approach  in  the  context  of  the  spectral 
multiscale  method  introduced  by  Y.  Efendiev  and  J.  Galvis  [20]. 

The  spectral  multiscale  method  is  designed  to  address  the  highly  oscilla¬ 
tory  problem  (8)  in  the  case  when  the  ratio  between  the  maximum  and  the 
minimum  values  of  A^  is  large  (high  contrast  problem).  The  main  idea  is 
to  complement  the  standard  MsFEM  approximation  space  Wh  =  Span  {0f} 
(where  are  the  MsFEM  basis  functions,  see  Section  3  above)  with  eigen¬ 
functions  that  correspond  to  small  eigenvalues.  More  precisely,  consider  the 
local  eigenvalue  problem 

— div  in  oji,  (28) 

with  homogeneous  Neumann  boundary  conditions,  where  Ui  :=  supp(0f)  is 
the  support  of  the  MsFEM  highly  oscillatory  basis  functions  0f  and  where 
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A^.  Note  that,  by  def- 


the  weight  is  dehned  by  A^ 


inition,  the  hrst  eigenvalue  A£=o  is  exactly  zero,  and  the  eigenvector 
associated  to  that  eigenvalue  is  constant  in  ooi.  We  then  select  the  eigenvec¬ 
tors  associated  with  an  eigenvalue  below  a  certain  threshold  r  >  0, 
and  construct  the  so-called  spectral  multiscale  basis  functions 


Xh  ■=  01 

The  spectral  multiscale  method  consists  in  performing  a  Galerkin  approxi¬ 
mation  of  (8)  on  the  approximation  space  W  :=  Span(x|^).  Note  that,  since 
0000  is  a  constant  function,  the  function  Xii=o  is  proportional  to  (pj.  Hence 
the  space  W  contains  the  standard  MsFEM  approximation  space  Span(0|). 

Our  project  is  to  use  the  Reduced  Basis  approach  to  more  efficiently 
compute  the  eigenfunctions  needed  to  build  W. 
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Summary 

We  report  here  on  the  work  performed  during  the  first  year  (october  2009  - 
October  2010)  of  the  contract  FA  8655-10-C-4002  on  Multiscale  problems  in 
materials  science:  a  mathematical  approach  to  the  role  of  uncertainty. 

The  bottom  line  of  our  work  is  to  develop  affordable  numerical  methods  in 
the  context  of  stochastic  homogenization.  Many  partial  differential  equations 
of  materials  science  indeed  involve  highly  oscillatory  coefficients  and  small 
length-scales.  Homogenization  theory  is  concerned  with  the  derivation  of  av¬ 
eraged  equations  from  the  original  oscillatory  equations,  and  their  treatment 
by  adequate  numerical  approaches.  Stationary  ergodic  random  problems 
(and  the  associated  stochastic  homogenization  theory)  are  one  instance  for 
modelling  uncertainty  in  continuous  media.  The  theoretical  aspects  of  these 
problems  are  now  well-understood,  at  least  for  a  large  variety  of  situations. 
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On  the  other  hand,  the  numerical  aspects  have  received  less  attention  from 
the  mathematics  community.  Standard  methods  available  in  the  literature 
often  lead  to  very,  and  sometimes  prohibitively,  costly  computations. 

In  this  report,  we  hrst  review  an  approach  popular  in  particular  in  the 
computational  mechanics  community,  which  is  to  try  and  obtain  bounds  on 
the  homogenized  matrix,  rather  than  computing  it.  Only  computations  of 
moderate  difficulty  are  then  required.  However,  we  will  show  that,  not  un¬ 
expectedly,  this  method  has  strong  limitations. 

We  will  next  introduce  a  class  of  materials  of  signihcant  practical  rele¬ 
vance,  that  of  random  materials  where  the  amount  of  randomness  is  small. 
They  can  be  considered  as  stochastic  perturbations  of  deterministic  materials, 
in  a  sense  made  precise  below.  We  will  adapt  to  such  a  case  the  well-known 
Multiscale  Finite  Element  Method  (MsFEM),  and  design  a  method  which  is 
much  more  affordable  than,  and  as  accurate  as,  the  original  method. 

The  works  described  below  have  been  performed  by  Claude  Le  Bris  (PI), 
Frederic  Legoll  (Co-PI),  and  Florian  Thomines  (first  year  Ph.D.  student). 


1  Introduction 

Many  partial  differential  equations  of  materials  science  involve  highly  oscilla¬ 
tory  coefficients  and  small  length-scales.  Homogenization  theory  is  concerned 
with  the  derivation  of  averaged  equations  from  the  original  oscillatory  equa¬ 
tions,  and  their  treatment  by  adequate  numerical  approaches.  Stationary 
ergodic  random  problems  are  one  of  the  most  famous  instances  of  mathe¬ 
matical  uncertainty  of  continuous  media.  However,  the  elaborate  tools  and 
techniques  of  (i)  mathematical  probability,  stochastic  analysis,  and  (ii)  nu¬ 
merical  analysis  and  large-scale  computing  have  not  yet  permitted  practical 
computations.  These  are  most  often  accomplished  otherwise  by  the  engineer¬ 
ing  community,  using  more  traditional  approaches.  Despite  dehnite  achieve¬ 
ments  by  leading  experts,  numerical  analysis  of  stochastic,  and  more  gener¬ 
ally  speaking  non  periodic,  homogenization  problems  remains  in  its  infancy. 

The  purpose  of  this  report  is  to  present  the  recent  progress  we  have  made 
during  last  year  on  this  topic,  with  the  aim  to  make  numerical  random  ho¬ 
mogenization  more  practical.  Because  we  cannot  embrace  all  difficulties  at 
once,  the  case  under  consideration  here  is  a  simple,  linear,  scalar  second  or¬ 
der  elliptic  partial  differential  equation  in  divergence  form,  for  which  a  sound 
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theoretical  groundwork  exists.  We  focus  here  on  different  practical  compu¬ 
tational  approaches. 

This  report  begins,  in  Section  2,  with  a  brief  introduction  to  stochastic 
homogenization.  There  is  of  course  no  novelty  in  such  an  introduction,  the 
only  purpose  of  which  is  the  consistency  of  this  report  and  the  convenience  of 
the  reader  not  familiar  with  the  theory.  We  will  recall  there  why  stochastic 
homogenization  often  leads  to  extremely  expensive  computations. 

In  Section  3,  we  describe  a  classical  approach  from  the  applied  commu¬ 
nities,  which  is  to  try  and  obtain  bounds  on  the  homogenized  matrix,  rather 
than  computing  it.  The  computational  gain  is  evident.  We  will  report  on 
some  numerical  experiments.  Such  experiments  are  likely  to  not  be  new.  But 
they  at  least  show,  quantitatively  and  qualitatively,  the  strong  limitations  of 
such  an  approach. 

As  pointed  out  above,  random  homogenization  for  general  stochastic  ma¬ 
terials  is  very  costly.  Yet,  it  turns  out  that  it  is  possible  to  identify  classes  of 
materials  of  signihcant  practical  relevance,  where  stochastic  homogenization 
theory  and  practice  can  be  reduced  to  more  affordable,  less  computationally 
demanding  problems.  These  materials  are  neither  periodic  (because  such  an 
oversimplifying  assumption  is  rarely  met  in  practice),  nor  fully  stochastic. 
They  can  be  considered  as  an  intermediate  case,  that  of  stochastic  perturba¬ 
tions  of  deterministic  (possibly  periodic)  materials.  The  case  when  the  tensor 
describing  the  properties  of  the  material  is  the  sum  of  a  periodic  term  and 
a  small  random  term  is  an  instance  of  such  an  approach.  In  Section  4,  we 
show  that  we  can  adapt  to  that  particular  setting  the  well-known  Multiscale 
Finite  Element  Method  (MsFEM),  which  is  designed  to  directly  address  the 
highly  oscillating  elliptic  problem,  rather  than  studying  the  limit  problem 
when  the  typical  small  lengthscale  goes  to  0.  This  method  has  been  initially 
proposed  for  deterministic  problems  [24,  21,  22,  14],  and  has  been  recently 
adapted  to  the  stochastic  setting  [18].  It  then  leads  to  extremely  intensive 
computations.  We  show  in  the  sequel  that,  if  the  problem  is  only  weakly 
stochastic,  then  it  is  possible  to  design  a  method  as  accurate  as  the  original 
MsFEM,  with  a  much  smaller  computational  cost.  As  we  explain  below,  this 
method  is  accurate  provided  the  stochastic  perturbation  is  indeed  small. 

We  collect  in  Section  5  some  conclusions  about  the  work  performed  so 
far,  and  future  directions  for  the  next  two  years  of  contract. 


3 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


2  Basics  of  stochastic  homogenization 

[Detailed  presentation  can  he  read  in  [1].] 

Stochastic  homogenization  is  best  understood  in  the  light  of  the  easiest 
context  of  homogenization:  periodic  homogenization.  This  is  the  reason  why 
we  begin  with  Section  2.1  laying  some  groundwork  in  the  periodic  context, 
before  turning  to  stochastic  homogenization  per  se  in  Section  2.2. 

We  refer  to,  e.g.,  the  monographs  [15,  19,  25]  for  more  details  on  homoge¬ 
nization  theory,  and  to  the  review  article  [1]  that  we  wrote,  addressing  some 
computational  challenges  in  numerical  stochastic  homogenization.  A  super 
elementary  introduction  is  contained  in  [13]. 

In  this  section,  we  present  classical  results  of  the  literature.  The  reader 
familiar  with  stochastic  homogenization  can  proceed  directly  to  our  contri¬ 
butions,  detailed  in  Sections  3  and  4. 


2.1  Periodic  homogenization 

For  consistency,  we  recall  here  some  basic  ingredients  of  elliptic  homogeniza¬ 
tion  theory  in  the  periodic  setting.  We  consider,  in  a  regular  bounded  domain 
D  in  W^,  the  problem 


-div 


A. 


per 


X 


=  /  in  T>, 


=  0  on  dV, 


(1) 


where  the  matrix  Aper  is  symmetric  definite  positive  and  Z'^-periodic.  We 
manipulate  for  simplicity  symmetric  matrices,  but  the  discussion  carries  over 
to  non  symmetric  matrices  up  to  slight  modihcations. 

The  microscopic  problem  associated  to  (1),  called  the  corrector  problem 
in  the  terminology  of  homogenization  theory,  reads,  for  p  hxed  in 

I  -div(Aper(|/)  (p-h  Vwp))  =  0  in  W^, 

1  Wp  is  ZAperiodic. 

It  has  a  unique  solution  up  to  the  addition  of  a  constant.  Then,  the  homog¬ 
enized  coefficients  read 


[^1p 


+  Ap^r{.y)ejdy, 


(3) 
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where  Q  is  the  unit  cube,  and  where  Wei  denotes  the  solution  to  (2)  for  p  =  e^, 
with  Ci  the  canonical  vectors  of  The  main  result  of  periodic  homogeniza¬ 
tion  theory  is  that,  as  e  goes  to  zero,  the  solution  to  (1)  converges  to  u* 
solution  to 


— div  [A*Vu*]  =  /  in  "D, 
u*  =  0  on  dV. 


(4) 


Several  other  convergences  on  various  products  involving  Ap^r  y—j  and 
also  hold.  All  this  is  well  documented. 


The  practical  interest  of  the  approach  is  evident.  No  small  scale  e  is 
present  in  the  homogenized  problem  (4).  At  the  price  of  only  computing  d 
periodic  problems  (2)  (as  many  problems  as  dimensions  in  the  ambient  space) 
the  solution  to  problem  (1)  can  be  efficiently  approached  for  e  small.  A 
direct  attack  of  problem  (1)  would  require  taking  a  meshsize  smaller  than  e. 
The  difficulty  has  been  circumvented.  Of  course,  many  improvements  and 
alternatives  exist  in  the  literature. 


2.2  Stochastic  homogenization 

The  mathematical  setting  of  stochastic  homogenization  is  more  involved  than 
that  of  the  periodic  case. 

We  put  ourselves  in  the  usual  probability  theoretic  setting  for  stationary 
ergodic  homogenization,  with  the  exception  that  our  notion  of  stationarity 
is  discrete.  It  intuitively  means  the  following.  Pick  two  points  x  and  y  ^ 
X  at  the  microscale  in  the  material  and  assume  y  =  x  +  k  with  k  G 
The  particular  local  environment  seen  from  x  (that  is,  the  microstructure 
present  at  x)  is  generically  different  from  what  is  seen  from  y  (that  is,  the 
microstructure  present  at  y).  However,  the  average  local  environment  in  x 
is  assumed  to  be  identical  to  that  in  y  (considering  the  various  realizations 
of  the  random  material).  In  mathematical  terms,  the  law  of  microstructures 
is  the  same.  This  is  stationarity.  On  the  other  hand,  ergodicity  means  that 
considering  all  the  points  in  the  material  amounts  to  hxing  a  point  x  in  this 
material  and  considering  all  the  possible  microstructures  present  there. 

2.2.1  Main  result 

With  the  same  setting  as  that  described  for  periodic  homogenization,  we 
may  now  briefly  describe  the  main  result  of  stochastic  homogenization.  The 
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solution  to  the  boundary  value  problem 


-d.v(^(|.  c.)v«')=/  in  V, 

=  0  on  dV, 


converges,  when  e  — 0,  to  the  solution  u*  of  (4)  where  the  homogenized 
matrix  is  now 


[A*]ij  =  ^(^j  (e*  +  Vwei(l/,  ■)V  A  {y,  •)  Cj  dy^  . 

The  corrector  problem  now  reads 

— div  [A  {y,uj)  {p +  Vwp{y,u))]  =  0  on  W’’, 
Vwp  is  stationary,  ^  (^J  ')  =  0- 


(6) 


(7) 


A  striking  difference  between  the  stochastic  setting  and  the  periodic  setting 
can  be  observed  comparing  (2)  and  (7).  In  the  periodic  setting,  the  corrector 
problem  is  posed  on  a  bounded  domain,  namely  the  periodic  cell  Q.  In  sharp 
contrast,  the  corrector  problem  (7)  of  the  random  setting  is  posed  on  the 
whole  spaee  and  cannot  be  reduced  to  a  problem  posed  on  a  bonnded 


domain.  The  reason  is,  condition  ^  yj  ')  dy j  =  0  in  (7)  is  a  global 


lim 

R - >+oo  \Bjl 


Q 

condition.  It  indeed  eqnivalently  reads,  because  of  the  ergodic  theorem, 
1  f 

'  Vwp{y,  ■)  dy  =  0  for  any  seqnence  of  balls  Bji  of  radii  R. 

Br 

The  fact  that  the  random  corrector  problem  is  posed  on  the  entire  space 
has  far  reaching  consequences  for  numerical  practice.  Actually,  this  is  proba¬ 
bly  the  main  source  of  all  the  praetical  dijfieulties  of  stochastic  homogeniza¬ 
tion. 


2.2.2  The  direct  numerical  approach 

Practical  approximations  of  the  homogenized  problem  in  random  homoge¬ 
nization  are  not  easily  obtained,  owing  to  the  fact  that  the  corrector  prob¬ 
lem  (7)  is  set  on  the  entire  space.  In  practice,  trnncations  have  to  be  con¬ 
sidered,  and  the  actnal  homogenized  coefficients  are  only  obtained  in  the 
asymptotic  regime. 
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Let  us  now  be  more  explicit.  In  practice,  the  matrix  A*  is  approximated 
by  the  matrix 

i^Nhj  i^)  =  +  A{y,uj)  (ej  +  Vw^{y,uj))  dy,  (8) 

\Qn\  Jq^  V  / 

which  is  in  turn  obtained  by  solving  the  corrector  problem  on  a  truncated 
domain,  say  the  cube  Qn  =  {—N,  NY  C  W^: 

I  -div  (yl(-,o;)  (p  + Vwf(-,o;)))  =  0  on  R'^, 

1  wY{-,oj)  is  (^AT-periodic. 

Although  A*  itself  is  a  deterministic  object,  its  practical  approximation 
is  random.  It  is  only  in  the  limit  of  inhnitely  large  domains  Qat  that  the 
deterministic  value  is  attained  (the  convergence  limAr^oo  has 

been  shown  in  [16,  Theorem  1]). 

At  hxed  iV,  the  approximate  homogenized  matrix  is  random:  a  set  of 
M  independent  realizations  of  the  random  coefficient  A  are  therefore  consid¬ 
ered.  The  corresponding  truncated  problems  (9)  are  solved,  and  an  empirical 
mean  of  the  truncated  coefficients  (8)  is  inferred.  This  empirical  mean  only 
agrees  with  the  theoretical  mean  value  of  the  truncated  coefficient  within  a 
margin  of  error  which  is  given  by  the  central  limit  theorem  (in  terms  of  M). 
For  a  sufficiently  large  truncation  size  iV,  this  truncated  value  is  admittedly 
the  exact  value  of  the  coefficient.  The  overall  computation  described  above 
is  thus  very  expensive,  because  each  realization  requires  a  new  solution  to 
the  problem  (9)  of  presumably  large  a  size  since  N  is  taken  large. 


3  Bounds  for  homogenization 

[Work  expanded  in  [1].] 

Given  the  above  computational  workload,  practitioners,  especially  scien¬ 
tists  from  the  applied  communities  (computational  mechanics,  ...),  some¬ 
times  choose  to  avoid  computing  actual  homogenized  equations  and  concen¬ 
trate  on  bounds  on  the  homogenized  matrices  A*.  In  [1],  we  have  carefully 
studied  this  approach,  which  has  some  (strong,  as  will  be  seen  below)  limi¬ 
tations. 
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We  consider  here  the  specihc  case  of  composite  materials  consisting  of 
only  two  phases.  We  denote  by  A  and  B  the  associated  matrix  coefficients, 
modelling  the  properties  of  the  phases.  We  also  £x  the  average  volnme  frac¬ 
tion  0  of  the  phase  A.  For  simplicity,  we  assnme  here  that  0  is  nniform  over 
the  whole  material.  The  problem  is  to  find  all  possible  homogenized  materi¬ 
als,  that  is,  mathematically,  matrices  A*,  that  can  be  attained  homogenizing 
such  phases  A  and  B  with  the  volume  fraction  6. 

In  this  specific  case,  some  bounds  on  the  homogenized  coefficients  may 
be  established.  Here,  we  present  one  example  of  such  bounds  (actually  the 
most  famous  one).  The  case  we  consider  is  a  scalar  equation  of  the  type  (1) 
with  a  matrix  coefficient  A^{x)  that  needs  not  be  periodic,  nor  stationary 
ergodic,  and  that  reads 

=  x%x)A  +  (1  -  X%x))B 

where  x^{^)  is  the  caracteristic  function  of  phase  A.  Obtaining  estimates 
on  A*  without  being  in  position  to  explicitly  compute  A*  at  a  reasonable 
computational  price  is  the  whole  interest  of  the  approach  by  “bounds” . 


3.1  The  Hashin-Shtrikman  bounds 

Based  on  the  density  of  the  matrices  obtained  by  periodic  homogenization 
in  the  set  of  matrices  obtained  by  arbitrary  homogenization,  it  is  possible 
to  derive  the  following  Hashin-Shtrikman  bounds  on  A*.  In  the  sequel,  we 
assume  B  >  A. 

Under  the  above  assumptions,  any  homogenized  matrix  A*  satisfies  the 
upper  bound 


A*p  ■  p  <  Bp  ■  p  +  6  min  \2p  ■  rj  +  (B  —  A)  ^rj  ■  rj  (1  —  6)h{p)]  (10) 

r?eIR‘i  ^  ^ 


for  any  p  G  where  h{p)  is  defined  by 


h{p)  =  min 


Ip  -  /cp 


Bk  ■  k 

Similarly,  any  homogenized  matrix  A*  satisfies  the  lower  bound 

A*p  ■  p  >  Ap  ■  p  {1  —  6)  max  \2p  ■  p  —  (B  —  A)~^p  ■  p  —  6g{p)]  ,  (11) 

r?eIR‘i 
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where  g{r])  is  defined  by 


aiv) 


If]  ■ 

max  — — - 

k&zd,kjto  Ak  ■  k 


Furthermore,  the  upper  bound  can  always  be  attained:  for  any  p  G  there 
exists  a  function  y,  Z'^-periodic  and  that  generally  depends  on  p,  such  that 
for  the  matrix  A*  obtained  by  periodic  homogenization  of 

A{-)=x{-)A  +  {l-x{-))B, 
e  e  e 

the  inequality  (10)  becomes  an  equality  (see  e.g.  [30]).  Likewise,  the  lower 
bound  (11)  can  always  be  attained.  We  have  summarized  in  [1,  Section  2.3.2] 
a  proof  of  the  Hashin  Shtrikman  bounds. 


Remark  1  Besides  the  Hashin- Shtrikman  bounds,  many  other  estimates  have 
been  proposed,  such  as  the  dilute  approximation,  the  self-consistent  method  [31] 
and  the  Mori  Tanaka  methods  [28].  They  are  all  based  on  the  fact  that  the 
problem  of  a  single  inclusion  in  an  infinite  material  (Eshelby  problem)  is  ana¬ 
lytically  solvable  [23].  Similarly  to  the  Hashin- Shtrikman  bounds,  the  spatial 
distribution  of  the  phases  is  not  taken  into  account  in  these  other  bounds. 
The  accuracy  of  these  estimates  and  bounds  strongly  depends  on  the  contrast 
between  A  and  B  and  the  volume  fraction  9  as  shown  on  Figure  1  below. 


3.2  Numerical  illustration 


We  consider  a  two-phase  composite  with  A  and  B.  We  denote  by  a  the 
scalar  conductivity  of  A  (respectively  b  the  conductivity  of  B)  with  a  <  b. 
We  denote  by  d  the  dimension,  and  by  6  the  volume  fraction  of  A. 

We  consider  the  case  of  the  random  checkerboard,  for  which  the  exact 
homogenized  matrix  is  known:  A*  =  a*  Id  =  Id.  In  this  simple  case, 
the  different  bounds  and  estimates  have  an  analytical  form:  the  homogenized 
coefficient  a*  is  bounded  from  below  by  the  harmonic  mean  (often  called  the 
Reuss  bound)  and  from  above  by  the  arithmetic  mean  (often  called  the  Voigt 
bound): 


1 

e/a  +  {i-e)/b 


<  a*  <  9a  -\-  {1  —  9)b. 
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These  bounds  are  also  called  Wiener  Bounds  or  Paul  bounds.  In  this  case, 
the  Hashin-Shtrikman  bounds  detailed  above  read  (see  e.g.  [25,  page  188]): 


\  da  +  6{b  —  a)  J 


<  a*  <b 


d6{b  —  a) 
db  +  (1  —  6){a 


1 


and  the  Self-Consistent  model  leads  to  an  estimate  Agff  of  the  effective  con¬ 
ductivity  a*  defined  implicitly  (see  [26])  by 


Q  ^  '^eff  Q  _  0\  ^ 

a  2Aeff  b  -|-  2Aeff 


0. 


On  Figure  1  we  plot  these  bounds  and  estimates  for  different  values  of  the 
contrast,  defined  by  b/a,  for  a  =  1.  Note  that  in  this  case,  by  construction, 
the  volume  fraction  for  any  a  and  b  is  6  =  1/2.  In  Tab.  1,  we  collect  the 
values  of  all  these  bounds  and  estimates,  for  the  particular  case  a  =  1  and 
b  =  10. 


Figure  1:  Different  bounds  for  the  checkerboard  test  case. 

A*  Harmonic  HS-  SC  Model  HS-I-  Arithmetic 
3.16  m  2)38  TOO  4T9  TSO 

Table  1:  Values  of  bounds  and  estimate  for  a  contrast  of  &/a  =  10. 

We  verify  that,  for  the  critical  volume  fraction  6  =  0.5,  even  for  a  contrast 
which  is  not  extremely  large  (a  =  1  and  b  =  10),  the  range  of  homogenized 
matrix  atteignable,  given  by  the  Hashin-Shtrikman  bounds,  is  wide.  In  such 
a  case,  the  spatial  distribution  of  phases,  which  is  not  taken  into  account  on 
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the  bounds,  is  certainly  important.  Note  also  that  a  typical  case  for  real- 
world  composites  is  more  challenging  than  the  case  above,  since  the  contrast 
is  usually  larger  (of  the  order  100  rather  than  10)  and  the  volume  fraction  is 

similar. 

Our  numerical  example  therefore  shows  that,  in  many  cases,  the  Hashin 
Shtrikman  bounds  cannot  provide  accurate  estimates  of  the  homogenized  ma¬ 
trix.  For  a  contrast  of  10,  the  error  between  the  bounds  and  A*  is  larger  than 
25  %.  For  a  contrast  of  100,  the  upper  bound  is  three  times  as  large  as  the 
actual  homogenized  value,  which  is  itself  three  times  as  large  as  the  lower 
bound.  There  is  therefore  a  need  for  developing  efficient  numerical  methods 
that  provide  more  accurate  results. 

4  A  weakly-stochastic  MsFEM  approach 

[Work  expanded  in  [3].] 

In  this  section,  we  show  how  the  Multiscale  Finite  Element  Method  (Ms¬ 
FEM)  can  be  adapted  to  the  stochastic  setting.  We  refer  to  [20]  for  a  review 
on  the  MsFEM  approach.  Let  us  recall  here  that  this  method  is  designed 
to  directly  address  the  original  problem  (namely  (5)  in  the  case  of  interest 
here),  keeping  e  at  its  hxed  value,  rather  than  studying  the  limit  problem 
when  e  — 0  (as  we  do  in  Section  2,  going  from  (5)  to  (4)).  Another  interest 
of  this  method  is  that  it  does  not  require  any  explicit  formula  for  the  homog¬ 
enized  tensor  (such  as  (2)-(3),  or  (6)-(7)),  which  are  not  always  available. 
More  details  and  comprehensive  numerical  tests  are  published  in  [3].  See 
also  [4]. 

4.1  MsFEM  approaches 

For  consistency  and  to  set  our  notation,  we  briefly  review  the  classical,  de¬ 
terministic  setting  for  MsFEM  approaches.  We  next  turn  to  the  stochastic 
setting.  We  consider  problem  (1),  which  we  reproduce  here  for  convenience, 

f  — div(A^(a;)Vu^(a;))  =  f{x)  in  V,  ,  . 

(  =  0  on  dV, 

where  is  a  symmetric  matrix  satisfying  the  standard  coercivity  and  bound¬ 
edness  conditions.  Note  that  the  approach  is  not  restricted  to  the  periodic 
setting,  so  we  do  not  assume  that  A‘^{x)  =  A{x/e)  for  a  periodic  matrix  A. 
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As  recalled  above,  we  wish  here  to  keep  e  hxed  at  a  (small)  given  value. 
The  MsFEM  approach  consists  in  performing  a  variational  approximation  of 
(12)  where  the  basis  functions  are  dehned  numerically  and  encode  the  fast 
oscillations  present  in  (12).  In  the  sequel  we  will  argue  on  the  variational 
formulation  of  (12): 


Find  G  Hq{T>)  such  that,  Wv  E  Hq{'D),  Ae{u^,v)  =  b{v),  (13) 


where 


Ae{u,v)  = 


E 

*,1 


A 


du  dv 


(X 


IV 


dxi  dx. 


dx  and  h{v)  =  f  v  dx. 


IV 


We  introduce  a  classical  PI  discretization  of  the  domain  V,  with  L  nodes, 
and  denote  0°,  i  =  1, . . . ,  L,  the  basis  functions. 


Definition  of  the  MsFEM  basis  functions  Several  dehnitions  of  the 
basis  functions  have  been  proposed  in  the  literature  (see  e.g.  [24,  21,  22,  14]). 
They  give  rise  to  different  methods.  In  the  following,  we  present  one  of  these 
methods.  We  consider  the  problem 


/  -div(A"(a;)V0-’^) 

I 


=  0  in  K 
=  0°Ik  on(9K. 


(14) 


Note  the  similarity  between  (14)  and  the  corrector  problem  (2).  Note  also 
that  the  problems  (14),  indexed  by  K,  are  all  independent  from  one  another. 
They  can  hence  be  solved  in  parallel,  using  a  discretization  adapted  to  the 
small  scale  e. 


Macro  scale  problem  We  now  introduce  the  hnite  dimensional  space 

Wh  ■■=  span  {01,  i  =  1, . . .  ,L}  , 

where  0f  is  such  that  0f|j^  =  0{’^  for  all  K,  and  proceed  with  a  standard 
Galerkin  approximation  of  (13)  using  Wh'- 

Find  e  Wft  such  that,  'iv  E  Wh,  Ae{ul^,v)  =  b{v).  (15) 

The  dimension  of  Wh  is  equal  to  L:  the  formulation  (15)  hence  requires 
solving  a  linear  system  with  only  a  limited  number  of  degrees  of  freedom. 
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Numerical  illustration  In  order  to  illustrate  the  MsFEM  approach,  we 
solve  (12)  in  a  one  dimensional  setting  with 

A\x)  =  5  +  50sin^  , 

on  the  domain  V  =  (0, 1),  with  £  =  0.025  and  /  =  1000.  We  subdivide  the 
interval  (0, 1)  in  L  =  10  elements.  On  Figure  2,  we  plot  the  MsFEM  basis 
functions  in  a  reference  element  and  the  MsFEM  solution 


Figure  2:  Left:  Multiscale  basis  functions  0^’^  in  the  reference  element. 
Right:  MsFEM  solution  in  the  domain  (0, 1). 


Natural  adaptation  to  the  stochastic  setting  When  applied  to  the 
stochastic  problem 

f  — div(y4^(a;, a;)VM'^(a;,a;))  =  f{x)  in  "D,  .  . 

(  =  0  on  dV, 

where  the  practical  issue  is  to  build  an  estimate  of  the  mean  E(n^(a;,  •))  using 
a  Monte-Carlo  simulation  method,  the  natural  adaptation  of  the  MsFEM 
method  is  the  following:  for  each  random  realization  m,  hrst  construct  a 
MsFEM  basis  and  next  solve  the  macroscale  problem.  This  approach  requires 
a  signihcantly  large  number  of  computations,  since,  for  each  realization,  a 
new  basis  of  oscillating  functions  is  built,  and  a  problem  at  the  macroscale 
is  solved.  Such  an  approach  has  been  described  and  analyzed  in  e.g.  [18]. 

4.2  A  weakly  stochastic  setting 

As  seen  above,  considering  general  random  materials  lead  to  extremely  ex¬ 
pensive  computations.  The  question  arises  to  know  whether  this  general  ran¬ 
dom  context  is  really  relevant,  and  whether  adequate  modihcations  can  lead 
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to  substantial  improvements.  Our  line  of  thought  here  is  based  on  the  fol¬ 
lowing  two-fold  observation:  classical  random  homogenization  is  costly  but 
perhaps,  in  a  number  of  situations,  not  necessary.  Many  materials,  albeit 
not  deterministic,  are  not  totally  random.  Some  of  them  can  be  considered 
as  a  perturbation  of  a  deterministic  material.  The  homogenized  behaviour 
should  expectedly  be  close  to  that  of  the  underlying  deterministic  material 
(and  thus  tractable  from  the  practical  viewpoint),  up  to  an  error  depending 
on  the  amount  of  randomness  present. 

Model  We  introduce  and  study  here  a  specihc  model  for  a  randomly  per¬ 
turbed  deterministic  material  (we  refer  to  [12]  for  a  quick  overview  of  this 
setting,  and  some  of  the  associated  numerical  techniques,  and  to  [1]  for  a 
more  comprehensive  review  of  our  contributions  along  these  ideas).  We  are 
interested  in  the  following  elliptic  problem 

f  — div  (A^(a;,  cu)  Vu^)  =  f{x)  in  "D  C 
\  ufi  =  0  on  dV, 

that  is  (16)  with 

A^{x,u)  =  Al^{x,u)  =  Al{x)  +  riA\{x,u), 

where  77  G  M  is  a  small  parameter,  A^  is  a  deterministic  matrix  uniformly 
elliptic  and  bounded,  and  A\{x,uj)  is  a  bounded  matrix.  The  matrix  is 
hence  a  perturbation  of  the  deterministic  matrix  A^. 

Remark  2  The  above  setting  is  of  course  one  possible  setting  where  the  the¬ 
ory  may  be  developed.  Other  forms  of  random  perturbations  of  deterministic 
(possibly  periodic)  problems  could  also  be  addressed.  See  e.g.  [5,  6,  7,  8]  and 
the  review  article  [1]. 

In  the  case  (18),  a  MsFEM  method  alternative  to  the  one  presented  in 
Section  4.1  can  be  proposed.  The  idea  is  to  compute  the  MsFEM  basis 
functions  only  once,  with  the  deterministic  part  of  the  matrix  A^  and  next 
to  perform  Monte-Carlo  realizations  only  for  the  macro  scale  problems.  We 
refer  to  [3]  for  all  the  details. 

As  above,  we  hence  hrst  solve  (14),  with  A^  =  Ag,  and  build  the  hnite 
dimensional  space 


(17) 

(18) 


Wi  :=  span {</>’, 
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We  next  proceed  with  a  standard  Galerkin  approximation  of  (17)  using  Wft: 
for  each  m  G  {1, . . .  ,M},  we  consider  a  realization  and  compute 

E  yVh  such  that 

r  (9??  f 

WvEWh,  J  =  y  f{x)v{x)dx.  (19) 

Since  the  MsFEM  basis  functions  are  only  computed  once,  a  large  computa¬ 
tional  gain  is  expected,  and  this  is  indeed  the  case. 


Numerical  studies  We  now  estimate  the  performance  of  the  approach.  To 
this  aim,  we  compare  the  solution  of  the  above  approach  with  the  solution  of 
the  direct  approach  (of  Section  4.1)  and,  for  reference,  the  solution  to  (17) 
obtained  using  a  hnite  element  method  with  a  mesh  size  adapted  to  the  small 
scale  e.  Our  estimators  are  built  as  follows: 


e{ui,U2)  =  E 


I  |Mi  —  U2\\n 


N 


(20) 


where  N  is  the  norm  used,  Ui  and  U2  are  the  solutions  obtained  with  any  two 
different  methods.  The  expectation  is  in  turn  computed  using  a  Monte-Carlo 


method.  Considering  M  realizations  {X, 
\\ui{-,Uj)  -U2{-,Uj)\\n 


l<m<M 


of  a  random  variable 


(here  X{uj)  = 


||'W2(-,ci;)||v 

and  the  empirical  standard  deviation  (Tm  as 


we  compute  the  empirical  mean 


1 


M 


m=l 

M 


CO 


m=l 


As  a  classical  consequence  of  the  Central  Limit  Theorem,  it  is  commonly 
admitted  that  E(X)  satishes 


|E(X)-/iM(X)|<1.96 


We  consider  the  following  numerical  example.  Let  {Xk^i)  denote  a 

sequence  of  independent,  identically  distributed  scalar  random  variables  uni¬ 
formly  distributed  over  the  interval  [0, 1].  We  dehne  the  random  conductivity 
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matrix  as 


+  Psin(27ra;/e) 
+  Psin(27r|//e) 


2  +  sin(27r|//e)  \ 
2  +  Psin(27ra;/e)  / 


{1  +r]Xk,i{uj))ld2, 


with  the  parameters  P  =  1.8  and  e  =  0.025.  Then  we  compute  Uref  solu¬ 
tion  to  (17)  on  the  domain  V  =  (0, 1)^,  with  /  =  1.  Let  um  and  us  be 
its  approximation  by  the  general  MsFEM  approach  (of  Section  4.1)  and  the 
weakly-stochastic  MsFEM  approach  described  above,  respectively.  The  nu¬ 
merical  parameters  for  the  computation  are  determined  using  an  empirical 
study  of  convergence.  We  used  for  the  reference  solution  a  hne  mesh  of  size 
hf  =  e/40.  The  MsFEM  basis  functions  are  computed  in  each  element  K 
using  a  mesh  of  size  Hm  =  e/80.  The  coarse  mesh  size  is  H  =  1/30.  We 
consider  M  =  30  realizations. 

We  report  in  Tables  2  and  3  the  estimator  (20),  along  with  its  conhdence 
interval,  for  the  norms  and  L^(P),  respectively. 


V 

^ref^ 

e{us,UM) 

1 

0.1 

0.01 

8.12  ±0.19 

7.17  ±0.02 
7.15  ±0.002 

17.37  ±0.78 
7.62  ±0.07 
7.28  ±0.007 

15.51  ±0.87 

2.56  ±0.10 

1.39  ±0.002 

Table  2:  H^{'D)  error  (in 

%). 

V 

:  ^ref^ 

^ref^ 

e{us,UM) 

1 

0.1 

0.01 

0.56  ±0.08 

0.54  ±0.01 

0.53  ±0.001 

1.69  ±0.49 
0.57  ±0.06 

0.62  ±0.005 

1.47  ±0.50 
0.20  ±0.07 

0.11  ±0.003 

Table  3:  error  (in  %). 


We  observe  that,  when  t]  is  small  (here,  r]  <  0.1),  the  alternative  approach 
provides  a  function  us  that  is  an  approximation  of  Uref  as  accurate  as  um- 
Recall  that,  since  the  MsFEM  basis  has  only  been  computed  once,  the  cost 
for  obtaining  an  empirical  approximation  of  ]E(m5')  is  much  smaller  than  that 
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for  obtaining  the  corresponding  empirical  estimator  of  E(ma^).  This  demon¬ 
strates  the  efficiency  of  the  approach.  As  expected,  when  rj  is  not  small 
(say  ?7  ~  1),  the  accuracy  of  the  solution  us  computed  with  the  alternative 
approach  proposed  here  decreases. 


Elements  of  proof  In  [3],  we  have  analyzed  the  method  introduced  here 
in  the  one-dimensional  setting  (see  also  [4]).  For  the  sake  of  analysis,  we 

assume  that  the  highly  oscillating  coefficient  reads  a^{x,u)  =  ,a;j, 

where  satishes  the  standard  assumption  of  stochastic  homogenization  (see 
Section  2.2).  The  problem  (17)  now  reads 


(21) 


We  assume  that  the  randomness  is  small,  in  the  sense  (see  (18))  that 


ari{x,u)  =  Qperix)  +  7]  ai{x,Uj),  (22) 

where  a^er  is  a  deterministic,  periodic  function  and  rj  is  a  small  deterministic 
parameter. 

In  [3],  we  have  bounded  from  above  the  difference  between  u*,  the  so¬ 
lution  to  the  homogenized  equation  (4),  and  the  weakly-stochastic  MsFEM 
solution,  in  the  following  sense.  For  a  given  realization  of  the  random  co¬ 
efficient  ap{x,uj),  let  u{-,u)  be  the  weakly-stochastic  MsFEM  solution,  that 
solves  (19).  By  construction,  this  solution  is  a  linear  combination  of  the 
highly  oscillating  basis  functions: 

L 

u{x,uj)  =  ^Ui{uj)(j)l{x). 

i=l 

Let  n„_MsFEM(a^,  1^)  be  the  associated  representation  in  terms  of  standard  PI 
basis  functions: 

L 

Vw-MsFem{x,uj)  =  y^Ui{u)(tPi{x). 

i=\ 

We  have  the  following  estimate: 

||m*  -  nw-MsFEM(-,cc;)  11^1(0,1)  <  C  (h  +  ^  +  (23) 
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where  C  is  a  deterministic  constant,  independent  of  h,  e  and  77,  and  C{ri)  is  a 
deterministic  function,  bounded  when  77  — 0.  In  the  above  bound,  <Jh{uj)  is 
a  random  number,  independent  of  77,  that  depends  on  e,  h  and  the  random 
term  ai{x,u)  in  (22). 

Let  us  comment  on  (23).  Assume  that  77  =  0,  i.e.  the  problem  (21)  is  a 
periodic  problem.  Then  our  method  is  identical  to  the  standard  deterministic 
MsFEM  method,  and  we  recover  from  (23)  the  classical  bound  known  in  that 
case,  namely 

\\u*  —  t^MsFEMlInpO,!)  ^  C  (^h  +  . 

Assume  now  that  ai  is  deterministic.  Then  our  method  is  not  exactly  the 
MsFEM  method,  since  we  do  not  take  into  account  Oi  to  build  the  highly 
oscillating  basis  functions.  In  that  case,  turns  out  to  vanish,  and  we 

infer  from  (23)  that 

\\u*  —  u„_msFem||hi(o,i)  ^  C  (^h  +  —  +  'rf‘C{r]^  . 

We  hence  observe  that,  provided  the  term  neglected  to  build  the  basis  func¬ 
tions  is  small  (namely  77  1),  we  obtain  a  similar  accuracy  as  with  the 

standard  MsFEM  method. 

A  similar  conclusion  holds  in  the  general  case  (22).  Note  also  that  the 
bound  (23)  is  valid  for  any  realization  uj  of  the  randomness.  It  is  therefore 
a  more  precise  result  than  a  bound  on  the  expectation  of  the  error,  where 
all  random  realizations  are  averaged.  For  instance,  the  bound  (23)  allows  to 
understand  what  is  the  probability  distribution  of  the  error. 

5  Conclusions  and  plan  for  the  following  years 

In  this  report,  we  have  hrst  reviewed  an  approach  to  obtain  bounds  (here,  the 
Hashin-Shtrikman  bounds)  on  the  homogenized  matrix.  This  approach  only 
involves  computations  of  moderate  difficulty.  However,  we  have  outlined  the 
strong  limitations  of  such  a  strategy.  In  some  cases,  the  difference  between 
the  lower  and  upper  bounds  is  indeed  very  large.  The  obtained  estimations 
are  then  inaccurate.  This  motivates  the  development  of  efficient  numerical 
methods  that  provide  more  accurate  results. 

To  this  aim,  we  have  focused  on  weakly  stochastic  materials,  for  which  we 
successfully  adapted  the  well-known  Multiscale  Finite  Element  Method  (Ms¬ 
FEM).  We  have  proposed  a  method  with  a  much  smaller  computational  cost 
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than  the  original  MsFEM  in  the  stochastic  setting.  Provided  the  stochastic 
perturbation  is  indeed  small,  the  method  we  propose  is  as  accurate  as  the 
original  one. 

We  summarize  now  the  directions  of  research  we  wish  to  pursue  during 
the  following  years. 

A  variant  of  classical  random  homogenization  In  the  short  term,  our 
aim  is  to  study  a  particular  setting  for  stochastic  homogenization,  which  is 
not  the  classical  setting  described  in  Section  2.2  (where  the  random  coefficient 
A  in  (5)  is  stationary).  The  setting  we  wish  to  study  is  the  case  when  the 
random  coefficient  is  the  composition  of  a  standard  deterministic  and  periodic 
function  Aper  with  a  stochastic  diffeomorphism: 

A‘'(X,U)  =  Aper  (f’^) 

where,  for  any  random  realization  u,  the  application  x  h- >  <h(a;,ci;)  is  a  diffeo¬ 
morphism.  Formally,  such  a  setting  models  a  microstructure  that  is  periodic, 
in  a  given  reference  conhguration.  The  latter  is  only  known  up  to  a  certain 
randomness.  Materials  we  have  in  mind  are  ideally  periodic  materials,  where 
some  random  deformation  has  been  introduced,  for  instance  during  the  manu¬ 
facturing  process.  Othewise  stated,  these  are  periodic  materials  seen  through 
random  glasses!  This  setting  has  been  initially  introduced  in  [9],  where  the 
homogenized  problem  is  identihed. 

An  interesting  question  in  that  context  is  that  of  numerical  discretization. 
In  the  classical  context,  a  standard  procedure  is  to  solve  the  corrector  problem 
on  a  truncated  domain  (see  (8)  and  (9)).  The  convergence  of  the  procedure 
is  given  by  [16,  Theorem  1].  We  currently  work  on  a  similar  analysis  in  the 
context  of  (24)  (see  [2]). 

A  more  theoretical  question  is  to  precisely  understand  the  behaviour  of 
u^{x,uj)  —  u*{x)  when  e  goes  to  0,  where  is  the  solution  of  the  highly 
oscillating  problem  and  u*  the  solution  of  the  homogenized  problem.  In 
the  classical  setting,  and  in  the  one-dimensional  case,  the  convergence  of 
e~^^'^{u^{x,uj)  —  u*{x))  to  a  Gaussian  process  has  been  shown  in  [17].  This 
question,  in  the  context  of  (24),  is  addressed  in  [2]. 

The  setting  (24)  is  in  general  not  a  weakly  stochastic  setting,  as  the 
amount  of  randomness  present  in  <F  may  be  large.  Yet,  in  the  case  when  the 
diffeomorphism  $  is  close  to  the  identity,  namely 

^{x,u)  =  X  +  T]^{x,uj)  +  0{t]‘^)  (25) 
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for  a  small  deterministic  parameter  rj,  the  amount  of  randomness  turns  out 
to  be  small.  This  case  is  thus  another  instance  of  randomly  perturbed  deter¬ 
ministic  materials  (recall  Section  4.2,  where  we  introduced  another  weakly 
stochastic  setting,  and  Remark  2,  where  we  pointed  out  other  weakly  stochas¬ 
tic  settings).  The  case  (24)-(25)  has  been  studied  in  [10,  11]. 

A  Fast  Fourier  Transform  approach  In  the  course  of  our  investigations, 
we  have  identihed  the  following  tracks  of  research,  which  are  closely  related 
to  the  research  directions  of  the  contract. 

First,  in  the  periodic  homogenization  setting  recalled  in  Section  2.1,  a 
method  based  on  Fast  Fourier  Transform  has  been  proposed  in  [29,  27].  The 
idea  is  as  follows.  Let  Aq  be  a  constant  symmetric  positive  matrix.  The 
corrector  problem  (2)  is  equivalent  to 

f  -div(Ao(p  + Vwp))  =  div((Ape^(?/)  -  Ao)  (p-f  Vwp))  in 
1  Wp  is  Z'^-periodic. 

The  idea  of  [29,  27]  is  to  solve  this  problem  iteratively.  Knowing  the  iterate 
Wp  at  iteration  k,  the  next  iterate  Wp'^^  is  dehned  as  the  unique  solution  to 

I  -div  (Aq  (p  +  VwJ+^))  =  div  {{Aper{y)  -  Aq)  [p  +  VwJ))  in 
1  Wp'^^  is  Z^-periodic. 

As  Aq  is  a  tensor  independent  of  y  (in  contrast  to  Aper{y)),  the  above  problem 
can  be  solved  very  efficiently  using  a  Fourier  transform.  Hence,  rather  than 
solving  (2),  we  are  left  with  solving  many  times  a  simpler  problem. 

Our  aim  is  to  compare  this  iterative  method  with  the  standard  method, 
in  term  of  efficiency.  The  choice  of  Aq  is  most  probably  of  paramount  impor¬ 
tance,  since  the  convergence  rate  (and  also  the  fact  that  the  iterations  in  k 
converge  or  not)  depends  on  it.  We  have  already  run  some  preliminary  tests 
with  this  method,  but  dehnite  conclusions  are  yet  to  be  obtained. 

Second,  in  the  context  of  stochastic  homogenization,  approaches  using 
some  decomposition  of  the  random  matrix  A{x,  u)  in  (5)  would  be  worthwhile 
to  investigate. 
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